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Program  activities  tc  accomplish  these  objectives  were  organized  "’ider  five 
tasks.  Under  Task  1  -  Literature  Survey,  a  survey  was  performed  to  dev :rmine 
the  state-of-the-art  in  design  .and  analysis  of  bolted  composite  joints. 
Experimental  evaluations  of  joint  static  strength  were  performed  under  Tasks  2 
and  3.  In  Task  2  -  Evaluation  of  Joint  Design  Variables,  strength  data  were 
obtained  through  an  experimental  program  to  evaluate  the  effects  of  twelve 
Joint  design  variables.  In  Task  3  -  Evaluation  of  Manufacturing  and  Service 
Anomalies,  effects  of  seven  anomalies  on  joint  strength  were  evaluated  experi¬ 
mentally  and  compared  with  Task  2  strength  data.  Bolted  composite  joint 
durability  was  evaluated  under  Task  4  -  Evaluation  of  Critical  Joint  Design 
Variables  on  Fatigue  Life.  Seven  critical  design  variables  or  manufacturing 
anomalies  were  identified  based  on  Task  2  and  3  atten;; '1  data.  Under  Task  5  - 
Final  Analyses  and  Correlation,  required  data  reduction,  methodology  development 
and  correlation,  .and  necessary  documentation  were  performed. 

This  report  documents  all  program  activities  performed  under  Tasks  2,  3, 

4  and  5.  Activities  performed  under  Ta9k  1  -  Literature  Survey,  were  pre¬ 
viously  reported  on  AFFDL-TR-78-179. ^Static  strength  methodology  and  evaluations 
of  joint  static  and  fatigue  test  data  are  reported.  Analytic  studies  complement 
methodology  development  and  illustrate:  the  need  for  detailed  stress  analysis, 
i  the  utility  of  the  developed  "Bolted  Joint  Stress  Field  Model"  (BJSFM)  procedure, 
and  define  model  limitations.  For  static  strength  data,  correlations  with 
analytic  predictions  are  included.  Data  trends  in  all  cases  are  discussed 
relative  to  joint  strength  and  failure  node.  For  joint  fatigue  studies,  data 
trends  are  discussed  relative  to  life,  hole  elongation,  and  failure  mode 
behavior. 

This  final  report  is  organized  Ju  the  following  three  volumes: 

Volume  1  -  Methodology  Development  and  Data  Evaluation 

Volume  2  Test  Data,  Equipment  and  Procedures 

Volume  3  -  Bolted  Joint  Stress  Field  Model  (BJSFM)  Computer  Program 
User' 8  Manual 
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SECTION  I 


INTRODUCTION 


One  objective  of  this  program  was  to  develop  a  static 
strength  methodology  for  mechanically- fastened  composite  joints. 
Tliis  volume  documents  user-options  and  instructions  l'or  a  com¬ 
puter  program  to  analyze>  the  effects  of  stress  concentrations  on 
laminate  strength.  Entitled  "Bolted  Joint  Stress  Field  Model" 
(BJSFM',  it  computes  stress  distributions  on  a  lamina  or  laminate 
basis  for  unloaded  or  loaded  (bolt  bearing)  holes  in  isotropic  or 
anisotropic  materials.  Failure  predictions  based  on  lamina 
properties  and  one  of  several  failure  criteria  are  possible. 
This  volume  describes  the  formulation  and  input  data  requirements 
and  output  options.  Sample  problems  and  a  computer  program 
listing  are  also  included. 


SECTION  II 


PROGRAM  DESCRIPTION 


The  Bolted  Joint  Stress  Field  Model  has  been  developed  to 
facilitate  strength  analysis  of  isotropic  or  anisotropic 
materials  at  individual  fastener  holes.  Static  strength  of  an 
anisotropic  laminate  with  a  fastener  hole  is  predicted  using  a 
closed- form  analytic  approach  based  on  (1)  elastic  anisotropic 
theory  of  elasticity,  (2)  lamination  plate  theory  and  (3)  one  of 
several  optional  failure  hypotheses.  The  program  has  capability 
to  handle  strength  and  stiffness  anisotropy,  general  in-plane 
loadings,  as  shown  in  Figure  1,  multi-material  (hybrid)  laminates 
and  arbitrary  hole  (bolt)  sizes.  BJSFM  modular  substructuring  is 
illustrated  in  Figure  2.  Input  data  required  are:  lamina  mech¬ 
anical  properties,  in-plane  loadings,  hole  geometry,  and  hole 
loading.  Options  are  available  which  provide  computation  results 
after  each  program  block. 

The  stress  field  calculations  are  based  on  two-dimensional 
anisotropic  theory  of  elasticity  solutions  for  a  homogeneous, 
anisotropic  infinite  plate.  Laminate  stress  distributions  around 
an  unloaded  or  loaded  (bolt  bearing)  hole  are  calculated  using 
plane  stress  assumptions.  Laminate  stress  and  strain  distribu¬ 
tions  for  combined  bearing  and  bypass  loads  are  obtained  using 
the  principle  of  superposition  (Figure  3).  Fastener  bearing  in 
idealized  as  a  cosine  radial  stress  distribution  (Figure  4). 
Finite  width  corrections  for  loaded  holes  are  based  on  super¬ 
position  of  infinite  plate  results  as  shown  in  Figure  5. 
Infinite  plate  solutions  are  exact  while  corrections  for  finite 
width  joints  are  approximate  and  most  accurate  for  width- to- 
diameter  ratios  greater  than  four. 

Laminate  strains  are  calculated  using  material  compliance 
relations.  Laminate  compliance  coefficients  are  determined  using 
lamination  plate  theory  with  unidirectional  (lamina)  elastic 
constants,  lamina  orientations  and  thicknesses.  Strains  for 
individual  plies  along  lamina  principal  material  axes  are  calcu¬ 
lated  using  coordinate  transformations.  The  solution  is  strictly 
valid  only  for  homogeneous  media;  however,  it  has  been  assumed 
valid  for  mid-plane  symmetric  laminates. 

Laminate  failure  is  predicted  by  comparing  elastic  stress 
distributions  with  any  of  five  material  failure  criteria  on  a 
ply-by-ply  basis.  Failure  can  be  assessed  at  any  location  in  the 
field  of  the  plate. 
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Figur*  1.  Genera!  Load  Condition.  Analyzed  Using  GJSTm 
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Figur*  2.  Bolted  Joint  Stress  Field  Modal  -  BJSFM 

Computer  Program  Flow  Chart 
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Figure  4.  Assumed  Cosins  Solt-Losd  Distribution 
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SECTION  III 


USER’S  INSTRUCTIONS 


Conversational  interactive  procedures  are  used  for  specify¬ 
ing  input  data  for  operation  of  the  BJSFM  program  which  is  pro¬ 
grammed  to  accept  free  formatted  input  data.  A  user  may,  after 
becoming  familiar  with  the  input  procedures,  elect  to  delete 
input  instructions  and  receive  only  question  marks,  identifying 
all  required  input  data.  Using  the  various  output  options,  users 
may  receive  as  much  or  as  little  data  as  desired.  The  nature  of 
required  input  data  is  dependent  on  user-selected  output  data 
options;  the  BJSFM  program  automatically  adjusts  its  input  data 
requirements  to  accommodate  each  output  data  option. 

Up  to  eight  different  ply  orientations  and  three  different 
materials  may  be  input.  For  each  ply  orientation,  a  correspond¬ 
ing  thickness  must  be  specified  as  well  as  material  for  hybrid 
laminates.  Ply  thickness  may  be  either  actual  or  a  normalized 
thickness.  A  mid-plane  symmetric  stacking  sequence  is  assumed. 
Zero  degree  plies  are  oriented  parallel  to  the  X-axis.  Nomencla¬ 
ture  applicable  to  tl^e  BJSFM  is  summarized  in  Figure  6.  Positive 
angles  are  measured  Counterclockwise  from  the  X-axis.  All  input 
data  units  must  be  consistent. 

Unidirectional  lamina  material  stiffness  properties  arc 
required  input  data  for  each  different  material  specified. 
Unidirectional  lamina  strength  allowables  for  each  material  are 
required  only  if  failure  analysis  is  to  be  performed.  If  the 
maximum  strain  material  failure  criterion  is  used,  lamina  strain 
allowables  must  be  input;  otherwise,  input  lamina  allowables  are 
in  terms  of  stress. 

Any  set  of  in-plane  far-field  stresses  may  be  applied  to  an 
infinite  anisotropic  or  isotropic  plate  (Figure  1).  Bearing 
stress  direction  is  independent  of  far-field  stress  directions. 

The  BJSFM  is  only  capable  of  handling  finite  widths  for  bolt 
bearing  problems;  the  width,  W,  is  defined  as  perpendicular  to 
the  bolt  load  direction  (Figure  1).  The  stresses  calculated  in  a 
finite-width  bolt  bearing  problem  are  approximate  and  most 
accurate  for  width-to-diameter  ratios  greater  than  four.  In 
combined  loading  conditions,  the  finite  width  routine  applies  to 
only  the  loaded  hole  portion  of  the  problem.  To  obtain  infinite 
plate  results  for  a  loaded  hole,  input  specimen  width  as  0.0. 

The  user  must  specify  the  "range"  (between  low  and  high)  of 
angular  interval  between  locations  around  the  hole  for  which  data 
will  be  calculated.  This  range  must  also  be  su, divided  by  user- 
selected  "degrees  between  output"  to  specify  points  at  which 
calculations  are  to  be  made.  "Step  increments"  are  used  to 
obtain  data  at  increasing  distances  away  from  the  hole  boundary. 
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Coordinate  System* 

X  -  V  Laminate  Axis  System  Originating  at  the  Center  of  the  Hole.  Zero  Degree  Plies  are  Parallel  to  the 
X-Axis. 

X'  -  Y'  Rotation  of  the  X-Y  Axis  System  for  Application  of  Far-Field  Stresses 
1  -  2  Lamina  Axis  System;  Fibers  are  Parallel  to  the  1  Axis  and  Transverse  to  the  2  Axis. 

Variable  Description  (Units)** 

2 

El  Lamina  Modulus  of  Elasticity  in  Fiber  (1)  Direction  (F/L  ) 

o 

E2  Lamina  Modulus  of  Elasticity  in  Transverse  (2)  Direction  (F/L  ) 

G12  Lamina  Shear  Modulus  (F/L2) 

VI 2  Lamina  Poisson's  Ratio 

EX  Laminate  Modulus  of  Elasticity  in  X  Direction  (F/L  ) 

EY  Laminate  Modulus  of  Elasticity  in  Y  Direction  (F/L2) 

GXY  Laminate  Shear  Modulus  (F/L2) 

VXY  Laminate  Poisson's  Ratio 

o 

T1  Lamina  Allowable  Tensile  Strength  in  Fiber  (1)  Direction  (F/L  or  L/L) 

Cl  Lamina  Allowable  Compressive  Strength  in  Fiber  (1)  Direction  (F/L2  or  L/L) 

T2  Lamina  Allowable  Tensile  Strength  in  Transverse  (2)  Direction  (F/L2  or  L/L) 

C2  Lamina  Allowable  Compressive  Strength  in  Transverse  (2)  Direction  (F/L2  or  L/L) 

S  l  amina  Allowable  Shear  Strength  (F/L2  or  L/L) 

PX  Stress  in  X'  Direction  (F/L2)  -  Independent  of  Input  Thickness 

PY  Stress  in  Y*  Direction  (F/L2)  -  Independent  of  Input  Thickness 

PXY  Shear  Stress  (F/L2)  -  Independent  of  Input  Thickness 

P  Applied  Bearing  Stress  (F/L2)  •  P  =  Bolt  Load/(Dia  x  Actual  Thickness) 

U  Displacement  in  X  Direction  (Lj 

V  Displacement  in  Y  Direction  (L) 

W  Specimen  Width  (L)  -  Bolt  Loading  Only 

DIST  Radial  Distance  from  Hole  Boundary  (L) 

a  Angle  of  Applied  Bolt  Load  with  X  Axis  (i.e.,  Bolt  Loading  Angle) 

(3  Rotation  Angle  of  X'  -  Y'  Axes  from  X  -  Y  Axis  System  (i.e.,  Off  Axis  Loading  Angle) 

0  Angle  from  X  Axis  to  a  Point  Around  Fastener  Hole 

0  Rotation  Angle  of  1  -  2  Axes  from  X  -  Y  Axis  System  (i.e..  Ply  Orientation  Angle) 

All  angular  measurements  are  positive  counterclockwise  from  the  X  axis. 

•  See  also  Figure  1 . 

*•  Any  consistent  set  of  units  may  be  used: 

F  *  Force 
L  ■  Length 

OPIMllsat 

Figure  6.  Program  Nomenclature* 
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The  maximum  input  step  increment  io  seven  evenly  spaced  concen¬ 
tric  circles;  the  first  step  is  always  at  the  hole  boundary. 


The  option  to  use  any  one  of  five  different  failure  criteria 
has  been  programmed  into  the  BJSFM.  Failure  analysis  is  applied 
on  a  ply-by-ply  basis;  therefore,  only  unidirectional  (lamina) 
allowables  are  required  input  data.  Only  the  maximum  strain  cri¬ 
terion  requires  the  allowables  to  be  input  as  strains;  all  others 
use  stress  allowables.  Equations  for  each  of  the  programmed 
failure  criteria  are  given  below.  When  the  right  hand  side  of 
any  of  the  equations  exceeds  unity,  failure  has  been  predicted 
for  the  ply.  Tension  or  compression  stress/strain  allowables 
used  in  each  criteria  are  selected  automatically,  depending  on 
the  sign  of  individual  stress  field  components  being  evaluated. 


Maximum  Strain 


1 


Maximum  Stress 


1 


Tsai-Hill 


Modified  Tsai-Wu 


Hoffman 


An  example  printout  of  the  "conversational" 
request  input  data  is  shown  in  Figure  7 . 


language  used  to 


RNH 

IDO  YOU  WANT  INSTRUCTIONS? 

?YES 

SELECT  DESIRED  OUTPUT  FROM  THE  FOLLOWING  CASES. 

1  CARPET  PLOT  DATA 

2  LAMINATE  PROPERTIES 

3  LAMINATE  STRESSES 

4  LAMINATE  STRAINS 

5  CIRCUMFERENTIAL  fc  RADIAL  STRESSES/ STRAINS 

6  DISPLACEMENTS 

7  STRAINS  PER  PLY 

8  STRESSES  PER  PLY 

9  FAILURE  CRITERIA  PER  PLY 

10  AUTOMATIC  SEARCH  FOR  FAILURE 
72,3.4,5,6,7.8,9,10 

INPUT  NUMBER  OF  DIFFERENT  PLIES  TO  BE  INPUT  (8  MAX)  AND 
NUMBER  OF  DIFFERENT  MATERIALS  (3  MAX) 

?4.  1 

INPUT  THE  UNIDIRECTIONAL  MATERIAL  PROPERTIES  FOR  EACH  MATERIAL 
IN  THE  FOLLOWING  ORDER:  El,  E2,  G12,  POISSONS  RATIO 
718.85E6. 1.9E6, .85E6, .3 

INPUT  THE  UNIDIRECTIONAL  ALLOWABLES  FOR  EACH  MATERIAL 
IN  THE  FOLLOWING  ORDER:  Tl,  Cl,  T2.  C2,  SHEAR 
7230000, 320000. 28200,32300, 17300 
INPUT  THE  ANGULAR  ORIENTATION  OF  EACH  PLY 

70. . 45. ,-45. ,90. 

INPUT  THE  THICKNESS  OF  EACH  PLY 

7.5. . 2. .2. .1 

INPUT:  FAR  FIELD  STRESSES  PX.PY.PXY,  OFF  AXIS  ANGLE,  BEARING  S'.  ’ESS 

AND  BOLT  LOADING  ANGLE. 

710000 . .  0 . . 2500 . . 45 . . 50000 ..10. 

I-ipUT  WIDTH  (0.0  FOR  INFINITE  PLATE) 

71.5 

INPUT  BOLT  DIAMETER.  DEGREES  BETWEEN  OUTPUT,  LOW  RANGE.  HIGH  RANGE, 
STEP  INCREMENT  AND  NUMBER  OF  STEPS  DESIRED  (7  MAX) 

7. 25. 30. . 0.. 360.,  .02,2 

INPUT  THE  NUMBER  WHICH  CORRESPONDS  TO  THE  FAILURE  CRITERIA 
YOU  WISH  TO  USE 

1  MAXIMUM  STRAI N 

2  MAXIMUM  STRESS 

3  TSAI-HILL 

4  MODIFIED  TSAI-WU 

5  HOFFMAN 
73 

OPU+llfrl* 


Figure  7.  Exampls  Input  Instructions 
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SECTION  IV 


OUTPUT  OPTIONS 


Various  output  options  are  available  for  user  selection. 
The  user  may  select  any  or  all  of  the  following  options  by  input¬ 
ting  the  appropriate  number(s): 

1.  Carpet  Plot  Data 

2 .  Laminate  Properties 

3.  Laminate  Stresses 

4.  Laminate  Strains 

5.  Circumferential  and  Radial  Stresses/Strains 

6„  Laminate  Displacements 

7.  Strains  per  Ply 

8.  Stresses  per  Ply 

9.  Failure  Criteria  per  Ply 

10.  Automatic  Search  for  Failure 

A  brief  description  of  each  output  option  follows. 

Option  #1  -  Carpet  Plot  Data  -  The  carpet  plot  data  routine 
will  automat  leal  1  y  vary  the  layup  of  a  user  i.iput  0,/+$/90<* 
laminate  family  and  calculate  any  one  or  ell  of  the  other  output 
options  (2  through  10)  for  each  layup.  Sixty-six  different 
layups  are  automatically  calculated  in  thie  routine?  therefore, 
large  amounts  of  data  will  be  generated  when  using  this  output 
option. 

Option  #2  -  Laminate  Properties  -  Laminate  stiffness  proper¬ 
ties  are  calculated  using  the  unidirectional  material  elastic 
constants,  ply  angular  orientations  and  ply  thicknesses.  These 
properties  are  calculated  with  respect  to  the  X-Y  axes  and  are 
the  same  as  would  be  obtained  using  conventional  lamination 
theory  approaches . 

Options  #3  and  #4  -  Laminate  Stresses/Strains  -  Laminate 

stress  and  strain  distributions  are  available  as  output  at  points 
around  the  perimeter  of  the  hole  and  at  other  user-specified 
concentric  circles  about  the  hole  boundary.  Principal  stresses 
and  strains  are  also  calculated.  All  output  is  referenced  to  the 
X-Y  axes.  Points  are  located  by  the  radial  distance  away  from 
the  hole  boundary  and  the  angular  orientation  from  the  X-axis 
(Figure  8) . 

Option  #5  -  Circumferential  and  Radial  Stresscs/Strains  - 

Circumferential  and  radial  laminate  stresses  and  strains  are 
calculated  by  a  coordinate  transformation.  Output  is  in  polar 
coordinates . 
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LAMINATE  STRESSES 


D15T 

ANGLE 

X  STRESS 

Y  STRESS 

0.000 

0.00 

-62603.22 

44412.49 

0.000 

30.00 

-40151.39 

-013.39 

0 . 000 

GO.  00 

-6019. 06 

-29542.53 

o .  ooo 

90.00 

50351 .00 

-11003.21 

0.000 

120.00 

401 27.47 

16033.31 

L- 

150.00 

9612.46 

20047.34 

— 

1  ^ 

SHEAR 
STRESS 
-9.44 
-34070. 37 
-19090.77 
-9.44 
77703.47 


MAX. 


MIN. 


PRINCIPAL  PRINCIPAL 


44417.50 
10054.04 
4557.17 
5035* 


-67093.73 

-59074.32 

-4<\; 


DIRECTION 
.01 


OP1M11UM 


Figure  8.  Laminate  Street  Distribution  Output  Data  Option 


Option  #6  -  Laminate  Displacements  -  Displacements  for  each 
point  are  output  as  U  and  V,  which  are  displacements  in  the  X  and 
Y  directions  respectively.  Due  to  limitations  in  the  derivation, 
displacements  for  the  loaded  hole  case  shall  be  considered  accu¬ 
rate  only  within  approximately  three  times  the  diameter  of  the 
fastener.  Unloaded  hole  displacements  are  exact  throughout  the 
plate. 

Options  #7  and  8  -  Strains/Stresses  per  Ply  -  Strains  and 
stresses  per  ply  are  calculated  and  output  in  the  lamina  (1-2) 
coordinate  system.  Eacli  ply  is  identified  along  with  the  loca¬ 
tion  of  the  point  around  the  hole  for  which  stresses/strains  ar 
calculated  (Figure  9). 
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STRAINS  PER  PLY 


DIST 

ANGLE 

PLY 

STRAIN  I 

STRAIN  2 

SHEAR  STRA 

0.000 

0.00 

0.00 

-.007128 

.011127 

— • 000004 

0.000 

0.00 

45.00 

.001990 

.002001 

.010254 

0.000 

0.00 

-45.00 

.002001 

.001990 

-.010?^ 

0.000 

0.00 

90 . 00 

.011127 

-.non.Ti 

0.000 

30.00 

0.00 

-.003509 

y' 

Figure  9.  Ply-byPly  Strain  Distribution  Output  Data  Option 


Option  #9  -  Failure  Criteria  per  Ply  -  Failure  criteria  per 
ply  option  applies  the  user-selected  material  failure  criterion 
(max.  stress,  Tsai-Hill,  etc.)  to  each  ply  using  the  input 
material  allowables.  The  "failure  number"  obtained  as  output 
data  indicates  the  value  calculated  by  the  failure  criterion 
using  the  stress  or  strain  components  at  a  point.  A  failure 
number  equal  to  or  greater  than  one  predicts  ply  failure.  The 
program  automatically  selects  tension  or  compression  allowables 
depending  on  the  sign  of  individual  stress/ strain  components 
being  evaluated.  "Failure  ratios"  are  output  which  indicate  the 
relative  magnitude  of  contributing  stress  components  to  the 
overall  failure  number.  Therefore,  fai^re  can  be  assessed  as  to 
which  stress  component  is  most  significant.  These  failure  ratios 
are  in  terms  of  the  lamina  (1-2)  coordinate  system. 
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Option  #10  -  Automatic  Search  for  Failure  -  The  automatic 
aearcn  For  FaTXui  e  routine  will  search  over  a  uaer-epecif ied 
range  at  each  amular  increment  for  the  most  critical  single 
point  as  calculated  by  the  material  failure  criterion.  Search 
for  failure  is  onb  done  at  the  first  step  increment  away  from 
the  hole  boundary.  Therefore,  if  a  search  for  failure  is  to  be 
performed  at  the  boundary  of  the  hole,  the  step  increment  must  be 
input  as  0.0.  The  program  will  automatically  ratio  the  input 
streusi  field  until  first  ply  failure  is  predicted.  Output  is  the 
in- plane  stresses  at  vhich  failure  is  predicted  along  with  the 
angular  orientation  oi'  the  predicted  failure  location.  Failure 
numbers  are  also  given  for  all  other  plies  at  the  critical  ply 
failure  initiation  angle.  Failure  ratios  are  also  output. 
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SECTION  V 


PROGRAM  LIMITATIONS 


The  following  are  the  limitations  of  the  BJSFM  program. 

o  Strictly  valid  for  homogeneous  anisotropic  flat  plates 
and  assumed  valid  for  mid-plane  symmetric  laminates. 

o  Displacements  inaccurate  at  poir.ts  more  than  three  times 
the  hole  diameter  away  from  the  hole  boundary  for  loaded 
hole  cases. 

o  Stress  fields  inaccurate  for  width  - Lo-diameter  ratios 
less  than  four. 

o  Maximum  of  eight  different  ply  angular  orientations 
( input) . 

o  Maximum  of  three  different  materials  for  hybrid  laminates 
( input) . 

o  Maximum  of  seven  steps  away  from  the  hole  (output) . 


The  following  equation  must  be  satisfied  to  obtain  valid  output. 
[(High  Range) -(Low  Range) ]/ (Degrees  Between  Output)  £72 


The  following  data  must  be  input  as  integers: 

Output  Option  Numbers 
Number  of  Different  Plies 
Number  of  Different  Materials 
Material  Number 
Number  of  Steps 
Failure  Criteria  Number 


IDO  YOU  MANY  INSTRUCT  I  OtIS  7 
7TE8 

SELECT  DESIRED  OUTPUT  7 ROM  THE  POLLOWINO  CASES • 
1  CARPET  PLOT  DATA 
LAMINATE  PROPERTIES 
LAMINATE  STRESSES 
LAMINATE  STRAINS 

CIRCUMFERENTIAL  A  RADIAL  STRESSES/ STRAINS 
DISPLACEMENTS 
STRAINS  PER  PLY 
STRESSES  PER  PLY 
FAILURr  CRITERIA  PER  PLY 
AUTOMATIC  SEARCH  POR  PM  f (URE 


1 

3 

4 

5 

s 

7 
E 
* 

10 

74.6 

INPUT  NUMBER  OP  DIFFERENT  PLIES  TO  DE  INPUT  (9  MAX)  AND 
NUMBER  OF  DIFFERENT  MATERIALS  (1  MAX) 

- ?4,1-| 

J Input  the  unidirectional  material  properties  for  sac  material 

IN  THE  FOLLOUIIIO  ORDER*  Cl,  B2.  (111,  POISSONS  RATIO 
*—  71S.8SC6.1.9C6,  .85R6,  .3 

INPUT  THE  ANGULAR  ORIENTATION  OF  EACH  PLY 
-70. .45. ,-45. ,90. 


— {  INPUT 

*--7.43,  . 


PLY 


THE  THICKNESS  OP  EACH 

25,. JS,. 00  — - 

INPUT*  PAR  FIELD  STRESSES  PX.PY.PXY, 


V-  42%  -  0°  Plies;.  50%  •  ±45°,  8%  90° 


OTP  AXIS  ANGLE,  REARING  STRESS 

AND  BOLT  LOADING  ANGLE. 

70,0,0,0,50000. ,20. 

INPUT  WIDTH  (0.0  POR  INFINITE  PLATE) 

71.5 

INPUT  BOLT  DIAMETER,  DEGREES  BETWEEN  OUTPUT,  LOW  RANOE,  HIGH  RANOE, 
STEP  INCREMENT  AND  HUMBER  OF  STEPS  DESIRED  (7  MAX) 

7. 25. 30. ,0., 90.. .02,3 


r  — 

LAMINATE 

strains  (Output  Option  No.  4) 

DIST 

ANGLE 

X  STRAIN 

Y  STRAIN 

SHEAR 

MAX. 

MIN. 

DXRECTIC 

STRAIN 

PRINCIPAL 

PRINCIPAL 

O.COO 

-  0.00 

-.006927 

.007109 

■OOOOOl 

.007109 

-.006927 

-.00 

0.000 

30.00  I 

-.004117 

.001918 

-.012140 

.005705 

-.008005 

18.20 

0.000 

Lr6o.oo 

.001704 

-.004583 

-.oioun 

.00*4099 

-.008979 

-33.87 

0.000 

'-90.00-' 

.005512 

-.006573 

.000001 

.005512 

-.006573 

.01 

.020 

0.00 

- • 005642 

.001 ,6< 

.000253 

.005767 

-.005643 

-1.27 

.020 

30.00 

-.003419 

.001412 

-.009993 

.004546 

-.006554 

38.21 

.020 

60.00 

.001920 

-.005407 

-.009304 

.003762 

-.007443 

-32.82 

-  .020 

90.00 

.003843 

-.004336 

.001000 

.003949 

-.004442 

12.34 

-  .040 

0.00 

-.004698 

.004788 

.000437 

.004793 

-.004703 

-2.63 

.040 

30.00 

-.002855 

.001130 

-.000420 

.003795 

-.005520 

38.34 

.040 

40.00 

.001898 

-.004792 

-.006542 

.003231 

-.005125 

-31.44 

.040 

90.00 

.002845 

-.003118 

.001857 

.002906 

-.001359 

15.96 

displacements  (Output  Option  No.  6) 


DIST 

ANGLE 

U 

V 

O.COO 

0.00 

.000794 

-.000000 

0.000 

30.00 

.000943 

.000342 

0.000 

60.00 

.000691 

.000535 

0.000 

90.00 

.000215 

.000294 

.020 

0.00 

.000655 

-.000093 

.020 

30.00 

.000794 

.000300 

.020 

60.00 

.000518 

.000471 

.020 

90.00 

.000059 

.000234 

.040 

0.00 

.000538 

-.000104 

.040 

30.00 

.000668 

.000268 

.040 

60.00 

.000378 

.000419 

.040 

90.00 

-.000045 

.000192 

DO  YOU  WISH  TO  CONTINUE? 
7NO 

STOP 


opivoni-no 

Figure  10.  Locd*d  Hol«  Cast 
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IDO  YOU  t/ANT  INSTRUCTIONS? 

7YES 

SELECT  DESIRED  OUTPUT  FROM  TNI  FOLLOWING  CASF.S. 
1  CARPET  PLOT  DATA 
LAMINATE  PROPERTIES 

laminate  stresses 


2 
I 

4 

5 

6 
7 
0 
9 

10 
73.5 

INPUT  NUMBER  OF  DIFFERENT  PLIES  TO  BE  INPUT  (8  MAX)  AND 
NUMBER  OF  DIFFERENT  MATERIALS  (3  MAX) 


LAMINATE  STRAINS 

CIRCUMFERENTIAL  L  RADIAL  STRESSES/ STRAINS 

DISPLACEMENTS 

STRAINS  PER  PLY 

STRESSES  PER  PLY 

FAILURE  CRITERIA  PF.n  PLY 

AUTOMATIC  SEARCH  FOR  FAILURE 


THE  UNI  DIRECTIONAL  MATERIAL  PROPERTIES  FOR  EACH  MATERIAL 
IN  THE  FOLLOWING  ORDERS  EI#  Z2 ,  012,  POISSONS  RATIO 

-7 10. E6 , I0.E6.3 .05E6.  .3  — - 

INPUT  THE  ANGULAR  ORIENTATION  or  EACH  PLY 
70. 

INPUT  THE  THICKNESS  Or  EACH  PLY 


Lfi 

L?i 


Isotropic  Properties 
(e1  =  E2,  G  =  E1 


2(1  +  v) 


INPUT l  FAR  FIELD  STCESSES  PX.PY.PXY,  OFF  AXIS  ANGLE,  BEARING  STRESS 
AND  BOLT  LOAD i NG  ANGLE. 

? 10000. ,0.,0.,0.,0..0. 

INPUT  BOLT  DIAMETER.  DEGREES  BETWEEN  OUTPUT,  LOW  RANGE,  HIGH  RANGE, 
STEP  INCRBIENT  AND  NUMBER  or  STEPS  DESIRED  (7  MAX) 
?-25,13.,0..90.,0.,l 


BIST 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 


laminate  stresses  (Output  Option  No.  3) 


ANGLE 

0.00 

15.00 

30.00 

45.00 

60.00 

75.00 

90.00 


X  STRESS  Y  STRESS 


-.00 

-490.49 

-.41 

5001.62 

15003.65 

25480.07 

29993.51 


-10000.00 

-6031.64 

-1.22 

5001.62 

5001.22 

1830.02 

-.00 


SHEAR 
STRESS 
0.00 
If  .0.53 
.70 
-50H.62 
-86.  2.36 
-68. 9. 72 
-.00 


MAX. 

MIH. 

DIRECTION 

PRINCIPAL 

PRINCIPAL 

-.00 

-10000.00 

0.00 

-.00 

-7322.13 

15.00 

-.00 

-1 .67 

30.00 

10003.25 

-.00 

45.00 

20004.87 

•  OC 

-30.00 

27318.88 

.00 

-15.00 

29993.51 

-.00 

-.00 

CIRCUMFERENTIAL  AND  RADIAL  STRESSES  6  STRAINS  (Output  Option  No.  5) 


DIST 

ANGLE 

THETA 

P-ADIAL 

SHEAR 

THETA 

RADIAL 

SHEAR 

STRESS 

STRESS 

STRESS 

STRAIN 

STRAIN 

3TRAIN 

0.000 

0.00 

-10000.00 

-.00 

0.00 

-.001000 

.000300 

0.000000 

0.000 

15.00 

-7322.13 

-.00 

-.00 

-.000732 

.000220 

.000206 

0-000 

30.00 

-1-62 

-.00 

.00 

-.000000 

.000000 

.000000 

o.ooo 

45.00 

10003.25 

-.00 

.00 

.001000 

-.000299 

- . 000000 

0.000 

60.00 

20004.87 

.00 

.00 

.002000 

-.000599 

. 000562 

o.ooo 

75.00 

27318.88 

.00 

.00 

.002731 

-.000019 

.000767 

0.000 

90.00 

29993.51 

-.00 

-.00 

.002999 

-•0009(50 

.000000 

DO  YOU 

WISH  TO 

CONTINUE? 

7  NO 


STOP 


GP1W11S-221 

Figure  11.  Isotropic  Unloaded  Hole 
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IDO  YOU  WANT  INSTRUCTIONS?  '  ' 

7YES  ^ 

SELECT  DESIRED  OUTPUT  PROM  THE  FOLLOWING  CASES.  f 

1  CARPET  PLOT  DATA  N 

2  LAMINATE  PROPERTIES  \l  Afe, 

1  LAMINATE  STRESSES  ^  ”  W 

4  LAMINATE  STRAINS  1  .  , 

5  CIRCUMFERENTIAL  «.  RADIAL  STRESSES/ STRAINS  \L 

6  DISPLACEMENTS  ^  “ 

7  STRAINS  PER  PLY  'I 

8  STRESSES  PER  PLY  ^  \  90°  i 

9  FAILURE  CRITERIA  PER  PLY  \  \  / 

10  AUTOMATIC  SEARCH  **OR  FAILURE  /  Wj'jT 

710  *  \ 

INPUT  NUMRER  OF  DIFFERENT  PLIES  TO  OE  INPUT  (8  MAX)  \ 

AND  NUMBER  OF  DIFFERENT  MATERIALS  (3  MAX) 

•  74.1-1 

ItffW  niE  UNIDIRECTIONAL  MATERIAL  PROPERTIES  TOR  EACH 
MATERIAL  IN  THE  FOLIX3UING  ORDER  i  El.  E2,  012,  POISSONS  RATIO 
71B.05E6, 1.9E6, . 85EG , .3 

INPUT  THE  UNIDIRECTIONAL  ALLOWABLES  FOR  EACH  MATERIAL 
IN  THE  FOLLOWING  ORDER l  Tl.  Cl,  T2 ,  Cl.  SHEAR 
?2.1E5,3.2E5.2.f32E4,  23E4.1.73E4 

INPUT  THE  ANGULAR  ORIENTATION  OF  EACH  PLY 

70. . 45. ,-45, ,90. 

INPUT  THE  THICKNESS  OF  EACH  PLY 

7.5. . 2. .2,  .1 

INPUT!  FAR  FIELD  STRESSES  PX.PY.FXY,  OFF  AXIS  ANGLE,  BEARING  STRESS 
AND  DOLT  LOADING  ANGLE. 

V 10000 . , 500 . , 2500 . , 30 . , 35000 . . 10 . - 

INPUT  WIDTH  (0.0  FOR  INFINITE  PLATE) 

70.0 

INPUT  BOLT  DIAMETER,  DEGREES  BETWEEN  OUTPUT,  LOW  RANGE,  HIOH  RANGE, 
STEP  INCREMENT  AND  NUMBER  Or  STEPS  DESIRED  (7  MAX) 

7 . 25 . 5 . .  0 . . 360 . . .02 . 2 

INPUT  THE  NUMBER  WHICH  CORRESPONDS  TO  THE  FAILURE  CRITERIA 
YOU  WISH  TO  USE 

1  MAXIMUM  STRAIN 

2  MAXIMUM  STRESS 

3  TSAI-HILL 

4  MODIFIED  TSAI-WU 

5  HOFFMAN 


•Input  Stresses  are  Ratioed 
Uniformly  Until  First  Ply 
Failure  is  Predicted 


AUTOMATIC  SEARCH  FOR  FAILURE  (Output  Option  No.  10) 
FAILURE  STRESSES 

PX  PY  PXY  P 

17243.40  862.17  4310.85  60351.88  - 


0  1ST 

ANGLE 

HLY 

5*  XI  LURE 

failure  ratios 

NUMBER 

2 

SHEAR 

020 

335,00 

0.00 

.355 

210 

.05  3 

020 

335.00 

45.00 

.999 

.811 

-.101 

.537 

020 

335.00 

-45.00 

.730 

-.400 

.650 

.620 

020 

i 

335.00 

i 

90.00  — 

— 1.000 

|  .606 

-.075 

:788J 

—  Indicates  Ply  Failure 

* — Failure  Initiation  Location  Angle 

User  Specified  Distance  Away  from  Hole 
at  Which  Failure  is  Predicted 


■Depicts  Relative  Magnitude  of 
Stress  Components  to  Failure  Number 


Note:  The  ''search''  is  only  done  within  the  user  specified  ''range'  with  an  accuracy  ot  1/2  the 
"degrees  between  output"  in  locating  failure  initiation. 


aptsans-its 


Figure  12.  General  Loading  Condition 
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PROGRAM  LISTING 


n  o 


C 

C  BOLTED  JOINT  STRESS  FIELD  MODEL  (BJSFM) 

C 

C  THIS  PROGRAM  COMPUTES  LAMINATE  STRESS  AND  STRAIN  DISTRIBUTIONS 

C  AROUND  LOADED  AND  UNLOADED  FASTENER  HOLES.  THE  PRINCIPLE  OF 

C  SUPERPOSITION  IS  USED  TO  OBTAIN  STRESS/ STRAIN  DISTRIBUTIONS 

C  FOR  A  GENERAL  APPLIED  LOADING.  OPTIONAL  FAILURE  ANALYSIS 

C  ROUTINES  ALLOWS  LAMINATE  STRENGTH  PREDICTIONS  USING  VARIOUS 

C  MATERIAL  FAILURE  CRITERIA  AND  HYPOTHESES.  PROGRAM  WAS 

C  DEVELOPED  BY  J .M.OGONOWSKI  OF  MCDONNELL  AIRCRAFT  CO., 

C  ST. LOUIS,  MISSOURI, 

C 

c 

c 

PROGRAM  BJSFM( INPUT=100,OUTPUT=150,TAPE105=INPUT, TAPE5=INPUT, 
ITAPE10a=OUTPUT, TAPES) 

C 

COMMON/ OWE/ El (3) , E2 { 3 ) , G1 2 ( 3 ) ,  VI 2 ( 3 ) 

COMMON/ TWO/ I OUT ( 1 5 ) , NUMPLY , NUMMAT , ANG ( 9 ) , PLYTHK ( 3 ) , MAT ID ( 3 ) 
COMflON/THREE/  IANG,  ILOW ,  IHIGH,  STPINK,  NUMSTP 
COMMON/ FOUR/ PX, PY, PXY, P, PW, ALPHA, BETA, DIA, CORRECT 
COMMON/ FIVE/ FXT( 3 ) , FXC(3) , FYT ( 3 ) , FYC ( 3 ) , FXY ( 3 ) , IFAIL 
COMMON/ S IX /AI (3,3) 

COM MON/ SEVEN/ S( 3, 3) 

COMMON/ EIGHT/ STRESS( 3 ,7,73), STRAIN( 3, 7 , 73 ) 

COMMON/N I NE/ STR1 ( 3 , 7 , 7  3 ) , STR2 ( 3 , 7 , 7  3 ) , STR1 2 ( 3 , 7 , 7  3 ) 

INTEGER  ANS, ANS2 , YES 
REAL  IANG, ILOW  t IHIGH 
C 

DATA  YES/ 'YES' /.NO/' NO1/ 


10  CONTINUE 

OUTPUT(6) *  IDO  YOU  WANT  INSTRUCTIONS?' 

READ (5, 20) ANS 
20  FORMAT (A3) 

C 

IF( ANS. EQ. YES)  WRITE(6,30) 

30  FORMAT ( '  SELECT  DESIRED  OUTPUT  FROM  THE  FOLLOWING  CASES.'/, 

1‘  I  CARPET  PLOT  DATA'/'  2  LAMINATE  PROPERTIES'/'  3  ', 

2 'LAMINATE  STRESSES'/'  4  LAMI NATE  STRAI NS ’ / ’  5  CIRCUMF' , 

3  '  ERENTIAL  S.  RADIAL  STRESSES/ STRAINS' /  '  S  UIS PLACEMENTS', '  /  , 

4'  7  STRAINS  PER  PLY'/'  0  STRESSES  PER  PLY'/. 

5 1  0  FAILURE  CRITERIA  PER  PLY*/'  10  AUTOMATIC  SEARCH  FOR 

6 'FAILURE' ) 

C 

DO  40  L=1 , 15 

IOUT(L)-0 

40  CONTINUE 

READ ( 5 ,  ) (IOUT(L) ,L=1 , 10) 

C 

IF(ANS.EQ. YES)  WRITE(5, 50) 

50  FOP.MAT{ '  INPUT  NUMBER  OF  DIFFERENT  PLIES  TO  BE  INPUT  (3  MAX)’, 
1*  AND'/'  NUMBER  OF  DIFFERENT  MATERIALS  (3  MAX)') 

READ{ 5 ,  ) NUMPLY, NUMMAT 

C 

IF( ANS. EQ. YES)  WRITE( 5, GO) 

60  FOHMAT( *  INPUT  THE  UNIDIRECTIONAL  MATERIAL  PROPERTIES  FOR  ’, 

1 ' EACH  MATERIAL'/'  IN  THE  FOLLOWING  ORDER:  El,  E2 ,  G12,  ', 

2 'POISSONS  RATIO') 

READ (  5  ,  )  (F.1(L)  ,E2(L)  .  G 1  ?  (  L )  ,  VI 2 (  L)  ,L=1,  NUMMAT) 

C 
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C 

TF( ANS.EQ. YES. AND. PUTOUT( IOUT, 9) .GE.l. )  WRITE(6,70) 

70  FORMAT ( *  INPUT  THE  UNIDIRECTIONAL  ALLOWABLES  FOR  EACH  * , 

1* MATERIAL*/’  IN  THE  FOLLOWING  ORDER i  T1 ,  Cl,  T2,  C2,  SHEAR') 
IF(PUTOUT( IOUT.9) .GE. 1. )  READ(5,  ) ( FXT ( J ) , FXC (  J )  , FYT (  J  )  , 

1FYC(J) ,rXV(J) , J-l.NUMMAT) 

C 

INPUT  ANGULAR  ORIENTATION,  THICKNESS  AND  MATERIAL 
IDENTIFICATION  PER  PLY 

IF( ANS.EQ. YES)  WRITE{6,80) 

80  FORMAT \ '  INPUT  THE  ANGULAR  ORIENTATION  OF  EACH  PLY') 

READ( 5 ,  ) (ANG(L),L-1,NUMPLY) 

IF( ANS.EQ. YES)  WRITE(6,90) 

90  FORMAT ( '  INPUT  THE  THICKNESS  OF  EACH  PLY* ) 

READ(5,  ) (PLYTHK(J),J-1,NUMPLY) 

100  CONTINUE 

I F ( NUMMAT . NE . 1 )  GO  TO  120 
DO  110  L-l.NUMPLY 
MATID(L)-1 
110  CONTINUE 
GO  TO  140 
120  CONTINUE 

IF( ANS.EQ. YES)  WRITE{6,130) 

130  FORMAT ( '  INPUT  THE  MATERIAL  NUMBER  FOR  EACH  PLY') 

READ( 5 ,  ) (MATID(L) ,L-1,NUMPLY) 

140  CONTINUE 

IF( PUTOUT( I OUT, 3 ) . EQ .0.0)  GO  TO  200 

IF* ANS.EQ . YES. AND. PUTOUT( IOUT, 3 ) .GE.l.)  WRITE(6,150) 

150  FORMAT ( '  INPUT;  FAR  FIELD  STRESSES  PX,PY,PXY,  OFF  AXIS  ANGLE,', 

1*  BEARING  STRESS'/'  AND  BOLT  LOADING  ANGLE.') 

IF(FUTOUT(IOUT,3) .GE.l.)  READ { 5 ,  )  PX , PY , PXY , BETA , P , ALPHA 

IF( ANS.EQ. YES. AND. P.NE. 0.0)  WRITE(6, 160) 

160  FORMAT( '  INPUT  WIDTH  (0.0  FOR  INFINITE  PLATE)*) 

IF(  PlirOUT(  IOUT,  3  )  .GE.  1 .  .AND.  P.NE.  0.0)  READ(5,  )  W 
C 

170  CONTINUE 

IF(ANS.EQ.YES. AND.PUTOUT( IOUT, 3) .GE.l.)  WRITE(6,190) 

180  FORMAT ( '  INPUT  BOLT  DIAMETER,  DEGREES  BETWEEN  OUTPUT,  LOW  ' , 

1' RANGE,  HIGH  RANGE,'/'  STEP  INCREMENT  AND  NUMBER  OF  STEPS', 

2'  DESIRED  (7  MAX)’) 

C 

IF(PUTOUT( IOUT, 3) .GE. 1. )  READ( 5 ,  )  DXA, IANG, I LOW, IHIGH, 

1STPINK, NUMSTP 
C 

IF( NUMSTP. GT.- 7)  NUMSTP-7 
RANGE*( IHIGH-ILOW) /IANG 

IF( RANGE. GT. 72)  OUTPUT(6)'  DEGREES  BETWEEN  OUTPUT  TO  SMALL;  ', 

1'  TRY  AGAIN' 

IF( RANGE. GT. 72 )  GO  TO  170 
C 

BI,-P*DIA 
PW-0 . 0 

IF( P.NE. 0.0. AND. W.NE. 0.0)  PW-BL/ ( 2 . 0*W) 

WD-W/DIA 

IF( WD. LT. 4.0. AND. W.NE. 0.0)  OUTPUT(6)'  CAUTION:  WI DTH-TO-DI AMETSR ' 
1'  RATIOS  LESS  THAN  4  GIVE  ERRONEOUS  RESULTS' 
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IF( ANS.  EQ.  YES.  AND.  PUTOUT(  I  OUT,  9)  .GE.  1  .  )  WRITE(6, 190) 

190  FORMAT ( '  INPUT  THE  NUMBER  WHICH  CORRESPONDS  TO  THE 

1' FAILURE  CRITERIA  '/'  YOU  WISH  TO  USE*/'  1  MAXIMUM  STRAIN’/. 
2'  2  MAXIMUM  STRESS'/'  3  TSAI-UILL'/'  4  MODIFIED  ', 

3 ' TSAI-WU * / '  5  HOFFMAN') 

IF(PUTOUT(IOUT,  9)  .GE.  1.  )  READ(  5  ,  )  IFML 

IF( PUTOUT( IOUT, 10) .EQ.2. )  NUMSTP=2 
20(t  CONTINUE 

CARPET  PLOT  ROUTINE 

IF( PUTOUT( IOUT, 1 ) .NE. 2 . )  GO  TO  240 
ANGl=AHG( l ) 

ANG2=ANG{ 2) 

ANG3=ANG( 3 ) 

At  JG4=ANG  ( 4 ) 
pLYTl=PLYTHK( 1 ) 

PLYT2=PLYTHK(2) 

PLYT3=PLYTHK(3) 

PLYT4=PLYTHK(4) 

DO  290  JKI=1 , 11 
FORFI V=— . 1 
DO  270  I JK— 1 ,11 
ANG(  1  )=ANG1 
ANG( 2)=AMG2 
ANG( 3 )=ANG3 
AUG( 4 ) =ANG4 
PLYTHK ( 1 )=PLYT1 
PLVTH:'(2)=PLYT2 
PLYTtIK  (  3  )  =PLYT3 
PLYTHK(4)=PLYT4 
NUMPLY-4 

PLYTHK (  1  )  =  1  .  1-JKI*  .  1 
FORFTV=FORFTV+ . I 
CHF.CK=FORFIV+PLYTHK(  1 ) 

IF ( CHECK. GT. 1.0)  GO  TO  280 
PLYTHK (  2  )=FORFIV/2  . 

PLYTHK ( 3 ) =FORFI V/ 2 . 

PLYTHK(4)=1 . 0-PLYTHK( 1 )-FORFIV 

10= PLYTHK ( 1 ) *100 . 4 

14  5= PLYTHK (2)*200.4 

I 90= PLYTHK (4)*100.4 

IF(IO.NE.O)  GO  TO  210 

NUMPLY-NUMPLY-1 

PLYTHK ( 1 )= PLYTHK ( 2) 

PLYTIIK(  2  )  =  PLYTHK(  3  ) 

PLYTHK ( 3 ) =PLYTHK ( 4 ) 

ANG( 1 )=ANG( 2 ) 

ANG(2)=ANG( 3) 

ANG( 3 )=ANG( 4 ) 

210  CONTINUE 

IF( 145 .NE . 0)  GO  TO  220 
NUMPLY=HUMPLY- 2 

I F ( NIJM PLY . EQ .  1 )  PLYTHK ( 1  )  =PLYTHK( 4  ) 

IF( NUMPLY.EQ. 1 )  AHG( 1 )=ANG(4 ) 

PLYTHK (2 )=PLYTHK(4) 

ANG(2)=AKG(4) 

220  CONTINUE 

IF( 190. EQ.O)  NUMPLY=NUMPLY-1 

IF( PUTOUT( IOUT, 1) .EQ.2. )  WRITE(6,230)  10,145,190 
230  FORMAT(//25X, 'LAYUP:  ',13,'/'  ,13,'/'  ,13) 
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BRANCHES  TO  OTHER  SUBROUTINES 

240  IF(PUTOUT(IOUT, 2) .EQ.O. )  GO  TO  260 
ALPH-ALPHA 
CALL  ABD(ALPH) 

CORRECT-1 .0 

DUMMY* PUT OUT ( I OUT ,98) 

250  CONTINUE 

IF( PUTOUT{ I OUT, 3) . EQ.O . )  GO  TO  260 
CALL  LAMSTR 

IF(PUTOUT(IOUT, 7) .EQ.O. )  GO  TO  260 
CALL  PLYSTR(IFAIL) 

IF( PUTOUT( I OUT, 9) .EQ.O. )  GO  TO  260 
CALL  FAILURE 

IF( PUTOUT( IOUT, 10) .EQ.2.)  DUMMY-PUTOUT( IOUT, 99) 
IF( CORRECT.LT. .999. OR. CORRECT. GT. 1 .001 )  GO  TO  250 


260  CONTINUE 

DUMMY-PUTOUT  (  I  OUT  ,.90 ) 

IF( PUTOUT ( IOUT, 1 ) .NE. 2 . )  GO  TO  290 

270  CONTINUE 
280  CONTINUE 
290  CONTINUE 

OUTPUT (6 ) 1  DO  YOU  WISH  TO  CONTINUE?’ 
READ( 5, 300)  ANS2 
300  FORMAT (A3) 

IF (ANS2.NE.NO)  GO  TO  10 


STOP 

END 

ft******************************************************** 

FUNCTION  PUTOUT(IOUT, IN) 

DIMENSION  IOUT( 15 ) 

PUTOUT-O . 

DO  10  J-1,15 

IF(IOUT(J) .GE. IN)  PUTOUT=l. 

10  CONTINUE 

DO  20  J-1,15 

IF( IOUT( J ) . EQ. IN)  PUTOUT*2. 

20  CONTINUE 

IF( IN. EQ.98)  DATA=0.0 
IF(  IN. EQ. 99 )  DATA-1.0 

IF(DATA.EQ. 1.0. AND.IN.LT. 10)  PUTOUT-O.5 

RETURN 

END 
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SUBROUTINE  ABD  (ALPHA) 

CALCULATES  LAMINATE  PROPERTIES  FROM  LAMINA  INPUT 
COMMON/ONE/ El ( 3 ) , E2 ( 3 ) , G12 ( 3 )  , V12 ( i ) 

COMMON/TWO/ IOUT( 15 ) ,NUMPLY,NUMMAT, ANG(8) ,PI  YTHK(8) , MATID(8) 
COMMON/ S IX/ AI( 3*3) 

COMMON/ S EVE N/S( 3, 3) 


DIMENSION  V21 ( 3 ) , DIV( 3) , Qll ( 3 ) , Q22 ( 3 ) ,012(3) ,Q66(3) ,U1(3) , 
1U2 ( 3 ) , U3 ( 3 ) , U4 ( 3 ) ,U5 ( 3 ) , QBAR(8 . 3 , 3 ) , ZZ( 16) , Z( 16) ,Q(3, 3) , 
2AA(  3 , 3)  ,  A(  3  f  3) 


REDUCED  STIFFNESSES  FOR  EACH  MATERIAL 

DO  10  M-l.NUMMAT 
V21(M)-E2(M)*V12(M)/E1(M) 

DIV  ( M)  — 1 . -V 12 (M) *V21(M) 

Q11(M)-E1(M) /DIV(M) 
Q22(M)-E2(M)/DIV(M) 
Q12(M)-V12(M)*E2(M)/DIV(M) 
Q66(M)-G12(M) 

10  CONTINUE 

INVARIANT  PROPERTIES 


DO  20  M-l.NUMMAT 

U1(M)-(3.*Q11(M)+3.*Q22(M)+2.*Q12(M)+4.*Q66(M) )/8. 
U2(M)»(Q11(M)-Q22(M) )/2. 

U3(M)«(Qll(M)+Q22(M)-2.*Q12(M)-4.*066(M) )/B. 
U4(M)-(Q11(M)+Q22(M)+6.*Q12(M)-4.*Q66(M) )/8. 
U5(M)-(Q11(M)+Q22(M)-2.*Q12(M)+4.*Q66(M) )/8. 

20  CONTINUE 

DO  30  1-1,3 
DO  30  J-l,  3 
A(I,J)=0. 

AA(I, J)=0. 

30  CONTINUE 

TRANSFORMED  REDUCED  STIFFNESSES  PER  PLY 
THICK-0. 

DO  40  L-l.NUMPLY 

DEG-ANG(L) *3.1415926535/180.0 

M— MATID(L) 

QBAR(L,1, l)=Ul(M)+U2(M)*COS(2.*DEG)+U3(M)*COS(4.*DEG) 
QBAR( L, 1 , 2 )=U4 (M)-U3(M) *COS( 4 . *DEG) 

QBAR(L, 2 , 2 )-Ul (M)-U2 (M) *COS ( 2 . *DEG 1+U 3 (M) *COS(4 . *DEG) 
QBAR(L, 1,3)=.5*U2(M)*SIN(2.*DEG)+U3(M)*SIN(4.*DEG) 
QBAR(L,  2, 3)-.5*U2(M)*SIN(2  .  *DEG  )-U3  (M)'*SIN(4.  *DEG ) 
QBAR(L, 3, 3)-U5(M)-U3(M) *COS( 4 . *DEG ) 

QBAR(L, 2, 1)-QBAR(L, 1, 2) 

QBAR(L, 3 , 1 )-QBAR(L, 1,3) 

QBAR(L, 3, 2)=QBAR(L, 2, 3) 

THICK=PLYTHK(L)+THICK 
ZZ(L+1 )— THICK 
40  CONTINUE 

Z(D  —  l.*THICK/2.0 
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CALCULATE  A  MATRIX 

DO  70  1-1,3 
DO  60  J-1,3 
DO  50  L-l.NUMPLY 
Z(L+1 )— Z( 1 )+ZZ( L+l ) 

ZA^Z(L+1)-Z(L) 

A(I, J)-A(I. J)+QBAR(L,I,J)*ZA 
50  CONTINUE 

MATRIX  Q  AND  Q0  ARE  DUMMY  MATRICIES  USED  IN  CALCULATIONS  INVOLVING 
THE  MANIPULATION  OF  OTHER  MATRICIES 

Q  (  I  •  J ) —A ( I ,  J ) /THICK 
60  CONTINUE 
70  CONTINUE 

COMPUTE  A/THICK  INVERSE  MATRIX 
I STEP- 1 

CALL  INVERSE{Q,AI) 

LAMINATE  MID-PLANE  PROPERTIES 

EX1-1 . 0/ AI (1,1) 

EY1-1 .  O/'AI  (  2 , 2  ) 

VXYl— — EX1*AI (1,2) 

GXYl— 1 4  0/AI (3,3) 

IF{ PUTOUTt I OUT ,  2 )  . EQ. 2 • )  WRITE(6,90)  SX1, EYl, GXYl. VXYl 
80  FORMAT(/25X, 'LAMINATE  PROPERTIES'/'  EX  -  ' , E9 . 3 . 2X, ' EY  -  ', 
1E9.3.2X, 'GXY  »  ' . E9. 3. 2X, ' VXY  -  '.F5.3) 

CALCULATE  MATERIAL  PROPERTIES  FOR  OFF-AXIS  BOLT  LOAD 

TRANSFORMED  REDUCED  STIFFNESSES  PER  PLY 

THICK-0. 

ALPHA-ALPHA* 3 . 141 5926 53 5/ 190 . 0 

DO  90  L-l.NUMPLY 

DEG-ANG( L ) *3 . 1415926535/180 . 0 

DEG-DEG- ALPHA 

M— MATID(L) 

QBAR(L.  I, 1)-U1 (M)+U2(M)*COS(2.*DEG)+U3(M)*COS(4.*DEG) 

QBAR( L,  1 , 2 )— U4 ( M) -U3 (M) *COS( 4 . *DEG ) 

QBAR(L, 2. 2)-Ul(M)-U2(M)*COS(2.*DEG)+U3(M)*COS(4.*DCG) 

QBAR( L, 1, 3)-.5*U2(M)*SIN(2.*DEG)+U3(M)*SIN(4.*DEG) 
QBAR(L.2,3)-.5*U2(M)*SIN(2.*DEG)-U3(M)*SIN(4.*DEG) 

QBAR( L,  3 , 3 )— U5 ( M)-U3 ( M) *COS( 4 • *DEG ) 

QBAR(L,  2. 1)-QBAR(L, 1, 2) 

QBAR(  L,  3 . 1 ) -QBAR( L, 1.3) 

QBAR( L.  3 , 2 )— Q3AR( L, 2.3) 

THICK— PLYTHK ( L ) +THICK 
ZZ( L+l )— THICK 
90  CONTINUE 

Z(l)— l.*THICK/2.0 
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CALCULATE  AA  :1ATRIX 

00  120  1-1,3 
DO  110  J-l , 3 
DO  100  L-l.NUMPLY 
'/ (  T.+ 1  )-7(  1  )  4.7. 7.  (1,4-1  ) 

ZA— Z( L+l )  —  Z (  L ) 

AA{ I, J)-AA( I, J)+QDAR(L, I, J)*ZA 

100  CONTINUE 

0(1, J)-AA(I, J)/THICX 

no  continue 

120  CONTINUE 

COM  PUTS  AA/TMICK  INVERSE  MATRIX 
I STEP- 4 

CALL  I  DIVERSE  (Q.S) 

OFF-AXIS  LA.' II MATE  PROPERTIES 

EX2-1 . 0/S ( 1 , 1) 

EY2-1 . 0/S ( 2 , 2) 

VXY2— EX2*S(  1,2) 

CXY2-1 . 0/S( 3 , 3) 


RETURN 

sun 

to************************************ 

Sin  ROUT  I  UK  INVERSE  (X,XI) 

CALCULATES  THE  INVERSE  OF  A  3X3  MATRIX 

DIMENSION  X( 3, 3) ,  XI (3, 3) 

COMMON  I STEP 


r>ET»(X(  1  .  1  )*X(  2)*X(  3, 3)  )  +  (X(  l  ,  2)  *X(  2,  3)*X(  3, 1 )  )  + 

1  (X (1.3) *X( 2,1) *X( 3, ?))-(*( 1,3) *X( 2,2) *X( 3,1))- 

2  (X(1,1)*X(2,3)*X(3,2) )-(X(l,2)*X(2,l)*X(3, 3) ) 

I F ( DOT . CQ .0.0)  GO  TO  10 

XI  ( l , 1 ) - ( X ( 2 . 2 ) *X ( 3 , 3 ) -X ( 2 , 3 ) *X ( 3 , 2 ) ) /OCT 
XI(l,2)-(X(2, 3)*X(3.1)-X(2,1)*X(3. 3)  J/DLT 
XI  l 1 , 3 )  =  ( X( 2 , 1 ) *X ( 3 , ? )-X( 2 , 2 ) *X ( 3 , 1 ) ) /DET 
XI  ( 2 , 2  )- ( X ( 1 . 1 ) *X ( 3 , 3 ) -X( l , 3 ) *X ( 3 , 1 ) ) / DET 
XI ( 2 , 3 ) = ( X ( 1 , 2 ) *X ( 3 , 1 ) -X { l , 1 ) *X ( 3 , 2 ) ) / DET 
XI  ( 3 , 3 )  —  ( X ( l , 1 ) *X( 2 , 2 )-X (1, 2 ) *X( 2 , 1 ) ) /DET 
XI  ( 2 , 1 ) -( X ( 3 , 2 ) *X ( 1 , 3 )-X ( 1 , 2 ) *X ( 3 , 3 ) ) /DCT 
XI  ( 3,  l)-(X(l, 2)*X(2, 3)-X(2,?)*X(l, 3) )/DET 
:CI(3,2)-(X(2,1  )*X(l,  3)-X(l.  1)*X(2,3)  )/DET 
GO  TO  30 

10  PRINT  20, I: TCP 

20  FOR’ SAT  (  1  SUM  ROUTINE  INVERSE  CALCULATES  A  SINGULAR  MATRIX 
1  ‘  AT  STEP  *13.' 

30  CONTINUE 
RETURN 
END 
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SUBROUTINE  LAMSTR 

CALCULATES  THE  LAMINATE  STRESSES  AND  STRAINS  DUE  TO  A 
QENERiVL  INPLANE  LOADING  WITH  A  DOLT  LOAD 

COMMON  /TWO/ 1  OUT  (IS).  MUMFLY ,  NUMt  lAT ,  ANG  (  9 )  ,  PLYTHK  ( 9  ) ,  MATID  (  9  ) 
OOMMON/THREE/ IANG, I LOW , IHIQH, STPINK, NUMSTP 
COMMON /FOUR/ PX.PY. PXY, P, PW. ALPHA, BETA, DIA, CORRECT 
COMMON/ SIX/AI(3,3) 

COMMON/ SEVE  N/ S ( 3 , 3 ) 

COMMON/ EIGHT/ STRESS  '  3,7,73), STRAIN ( 3,7,73) 

REAL  IANG,  I  LOW,  HUGH 

DIMENSION  STR( 3,7,73),U(7,73)#V(7,73) . UX ( 7 , 73 ) , VY( 7 , 73 ) 

PX-CORRECT*PX 
PY«CORREOT*PY 
PXY-CORRECT*  PXY 
P-CORRECT*P 
FW-CORRECV*PW 


PI-3. 14159:;  6:535 
NUMPT— (  ( IHIGH-ILOW) /IANG)  +  1 
DO  10  J-l, NUMSTP 
DO  10  K-l, NUMPT 
U  ( J  ,  K )  — O  .  0 
V(J,K)-0.0 
DO  10  1-1.3 
STRESS (I,J,K)-0.0 
STRAIN(I, J,K)-0.0 
10  CONTINUE 

CALCULATE  UNLOADED  HOLE  STRESSES 

Ii?(  PX  .  EQ.  0 . 0 )  GO  TO  20 
BETAO-BETA 

CALL  UNLODED ( PX , D I A , A I , BETAO , STRESS , U , V ) 
20  CONTINUE 

IP( PY.EQ.0.0)  GO  TO  40 
BETA90— BETA+90 . 0 

CALL  UNLODED (PY,D IA. AI, BETA90 , STR, UX, VY) 

DO  30  J-l ,  NUMSTP 

DO  30  K-l , NUMPT 

U(J,K)-U(J,K)-tUX(J,K) 

V(J,K)-V(J,K)4VY(J,K) 

DO  30  1-1,3 

STRESS ( I , J . K ) -STRESS ( I , J , K )+STR ( I ,  J , K ) 

30  CONTINUE 

40  CONTINUE 

IF( °XY . EQ . 0 . 0)  00  TO  70 
BETA4 5-BETA + 4 5 . 0 

CALL  UNLODED  (  PXY,  DIA,  AI ,  BETA45,  STR,  UX,  VY) 
DO  50  J-l, NUMSTP 
DO  50  K— 1 , NUMPT 
U( J , K)— U( J, K)+UX  t J, K) 

V ( J , K ) —V ( J , K ) +VY  t  J , K ) 

DO  50  1-1,3 

STRESS ( I , J, K)— STRESS ( I , J . K)+STR( I ,  J, K ) 

50  CONTINUE 
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B ETA45* 9 ETA-45 . 0 

pxyn—pxy 

CALL  UNLODKO(  PXYN, DIA, AI ,  BETA45,  STR.UX, VY) 
DO  GO  J-l, NUMSTP 
DO  GO  ’>1,NU'1PT 
U(J,K)-U(J,K)+UX(J,K) 

V(J,K)-V(J,E)+VY(J.K) 

DO  GO  1-1, 3 

STRESS ( I « J,  N ) -STRESS!  I .  J,  K)  «-3TR(  I ,  J,  K ) 

60  CONTINUE 
C 

70  CONTINUE 

CALCULATE  LOADER  HOLE  STRESSES 

IF(  P.  F.Q.0.0)  GO  TO  100 
AL  TIAO-ALPHA 
PB*»P 

CALL  LOADED ( PB, DIA, 0, ALPHAO, STR, UX, VY) 

DO  BO  J*1 , NUMSTP 
DO  30  K-l.NUMPT 
U(J,K)-U(J,K)+UX(J,K) 

V(J.K)*V(J.K)+VY(J.K) 

DO  BO  1-1,3 

STRESS ( I , J, K) "STRESS ( I ,  J  , K )+3TR(  I ,  J, K ) 

□0  CONTI  NUE 

ALPHAO* ALPHA 

CALL  UNLODE D( PW ,  OIA, AI , ALPHAO, STR, UX, VY) 

DO  90  J-l, NUMSTP 
DO  90  K*1  ,  tIUMPT 
U(J,K)-U{ J,K)+UX<J,K) 

V  (  J ,  K ) *V ( J  ,  K ) +VY (  J ,  K ) 

DO  90  1-1,3 

STRESS ( I , J, K ) “STRESS ( I , J , K )+5TU( I ,  J ,  K ) 

90  CONTINUE 
100  CONTINUE 

I F  (  PUTOUT  (  I  OUT ,  3)  .EQ. 2. )  WRITE(G.llO) 

110  FORMAT! ///23X , ’LAMINATE  STRESSES'//’  BIST 
1‘  Y  STRESS  SHEAR  MAX.  MIN. 

239X, ' STRESS  PRINCIPAL  PRINCIPAL') 


ANGLE  X  STRESS', 
DIRECTION' /, 
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CALCULATE  PRINCIPAL  STRESSES 

IF(PUTOUT( IOUT, 3) .NE. 2. )  00  TO  140 
1X5  130  J  J— 1 ,  NUMSTP 
DO  130  NN-1 , NUMPT 

PRINA—( STRESS  ( 1 ,  J  J , NN) -STRESS  ( 2 ,  JJ , NN) ) * ( STRESS ( 1 ,  JJ, NN)  - 
1  STRESS (  2  ,  JJ ,  Nf!}  )  /  4  . 

,>R INA-SQRT(  PRTNA-f  STRESS  {  3 , J J , NN )  *STRESS(  3  ,  J J .  NN )  ) 

FRINl—( STRESS ( l  «  JJ»tlN)+E  ?RESS(  2  JJ  ,  NN)  )/2  . +PRINA 
PRIN2-(STRESS(1, JJ,NN)+S  i'RESS(  2  ,  JJ,  NN)  )/2.-PRINA 
TSTS— STRESS ( 1 , J J , NN ) -STRESS ( 2 , J J , NN) 

DTRCT-O. 

IF(TSTS.NE. 0 . )  DIRCT-. 5*A1A«( 2 . * STRESS ( 3 , JJ, NN) /TSTS) 
DIRCT-180.*DIRCT/3. 1415926535 

IF( PUTOUT( IOUT, 3 ) .NE. 2  » )  GO  TO  140 
ANGLE* ( NN- 1 ) * I ANG+ I LOW 
DIST»(JJ-1  )‘STPINK 

WRITE( 6, 120)  DIST, ANGLE, STRFSS( 1 ,JJ,NN) , STRESS ( 2 , JJ , NN) , 
1STR£GS( 3, JJ,NN) , PRIN1, PKIN2, DIRCT 
120  FORMAT ( F6 .  3  ,  F8 . 2 , 5F 1 1 . 2 ,  FO  .  2  ; 

130  CONTINUE 
140  CONTINUE 

IF( PUTOUT( IOUT, 4) .E0.2. )  WRITE(6, 150) 

150  FORMAT( ///28X, *  LAMINATE  STRAINS*//'  DIST  ANGLE  X  STRAIN*, 
1*  Y  STRAIN  SHEAR  MAX.  MIN.  DIRECTION* / , 

239X,' STRAIN  PRINCIPAL  PRINCIPAL*) 

CALCULATE  LAMINATE  STRAINS 

DO  160  JJ-l, NUMSTP 
DO  160  NN-1 , NUMPT 
DO  160  KK-1,3 
DO  160  MM— 1 , 3 

STRAIN ( KK . J J . NN ) - AI ( KK , MM ) ‘STRESS ( MM , J J , NN ) +STRAI N ( KK , J  J , NN ) 
160  CONTINUE 

CALCULATE  PRINCIPAL  STRAINS 

IF(PUTOUT(IOUT,  4)  .NF..  2.  )  GO  TO  1  >0 
DO  180  JJ-1, NUMSTP 
DO  180  NN-1, NUMPT 

PR IMA— ( STRAIN ( 1 , JJ , NN) -STRAIN ( 2 , JJ , NN ) )*( STRAIN ( 1 , J J , NN) - 
1  STRAIN (  2 ,  JJ ,  Nil )  )  / 4  . 

PRI NA-SQRT  (  PR  I NA  l-STRAI  N  (  3  ,  J  J ,  Ntl )  *  .  2  5  *  STRAI N  (  3  ,  J  J ,  NM )  ) 

PRIN1— ( STRAIN ( 1 , J J , NN )+ STRAIN ( 2 , JJ , N  M ) )/2.+PRINA 
PRIN2— (STRAIN( 1 , JJ ,NN] +STRAIN( 2 , JJ , NN) )/2 .-PRINA 
TSTS— STRAI  N(  1 ,  JJ ,  NN )  -  ‘STRAIN  (  2  ,  J  J ,  NN) 

DIRCT— 0 . 

IF ( TSTS . NE . 0 . )  DIRCT- . 5* ATAN ( 2 . * STRAIN ( 3 , JJ, NN) /TSTS) 

DIRCT— 180 .‘DIRCT/ 3 . 1415926535 
DIST— (JJ-1 ) ‘STPINK 

ANGLE— ( NN-1 ) *1 ANG+ILOW 

WRITE(6, 170)  DIST. ANGLE, STRAIN( 1 , JJ , NN) , STRAIN(  2  ,  J J ,  NN)  , 
1STRAIN( 2 . JJ.NN) , PR INI , PRIN2 , DIRCT 
170  FORMAT ( F6.3,FB.2,5Fll.6,F8.2) 

180  CONTINUE 
190  CONTINUE 
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CALCULATE  CIRCUMFERENTIAL  AND  RADIAL  STRESSES  {.  STRAINS 


IF( PUTOUT! IOUT,  5) .EQ. 2. )  WR1TE(6, 200) 

200  FORMAT! ///15X, 'CIRCUMFERENTIAL  AND  RADIAL  STRESSES  (.  STRAINS', 

I//'  DIST  -ANGLE  THETA  RADIAL  SHEAR  THETA', 

2'  RADIAL  SHEAR’ /21X, ' STRESS  STRESS 

3 ' STRESS  STRAIN  STRAIN  STRAIN') 

IF! PUTOUT( IOUT, 5 ) .NE. 2. )  GO  TO  230 
DO  220  J-l.NUMSTP 
DO  220  N-l.NUMPT 

ENERGY- . 5* ( STRESS! 1 , J , N ) ‘STRAIN! 1 ,  J ,  N ) +STRESS! 2 , J ,  N ) * 

I  STRAIN! 2 , J, N )+ STRESS! 3 , J , N ) ‘STRAIN! 3 , J , N ) ) 

ANGLE- ! N- 1 ) * I ANG+ I LOW 
D-ANGLE*PI/ 180  0 
DIST— (J- I ) ‘STPINK 

RADSTS-STRESS ( 1 , J , N ) ‘COS ( D ) ‘COS ( D ) +STRESS (2,J,N)*SIN(D)‘ 

1  SIN!D)+2. ‘STRESS! 3,J,N)*SIN(D) ‘COS! D) 

CIRSTS-STRESSU,  J,N)*SIN!D)-SIN(D)  +  STRESS(2, J,N)‘COS!D)‘ 

1  COS! D)-2. ‘STRESS! 3, J, N ) *SIN( D ) *COS( D) 

SHRSTS— 1 . *  STRESS! I . J , N ) *S IN( D ) ‘COS ( D ) +STRESS ( 2 , J . N ) ‘SIN ( D ) * 

I  COS! DJ+STRESS! 3 , J , N ) * !COS( D) ‘COS (D)-SIN(D)*SINfD) ) 

RADSTN— STRAIN! I , J, N ) ‘COS !D) ‘COS( D )+STRAIN( 2,J,N)*SIN(D)‘ 

1  SIN ( D ) +STRAIN( 3 , J , N ) *SIN( D ) *COS( D ) 

Cl  RSTN— STRAIN!  li>J,N)*SIN!D)  ‘SIN!  D  )  +  STRAIN(  2  ,  J  ,  N  )  ‘COS  (  D )  * 

I  COS! D)-STRAIN( 3 , J ,N) *SIN( D) ‘COS! D) 

SHRSTN—1. ‘STRAIN!  1,J,N)«SIN!D)  *COS(  D)+STRAIN(  2,  J,N)‘SIN(D)‘ 

1  COS!D)+STRAIN( 3. J,N)*(COS(D)*COS(D)-SIN(D)‘SIN(D) ) 

IfiUTE! 6,210)  DIST, ANGLE , Cl RSTS , RADSTS , SHRSTS , Cl RSTN , RADSTN , 

1  SHRSTN 

210  FORMAT! F6.3,F10.2,F12.2,2F11.2, 3F9 . & ) 

220  CONTINUE 
230  CONTINUE 

OUTPUT  DISPLACEMENTS 

U3— 3 . 0*DIA 

DISP-DIA/ 2 • 0+NUMSTP*STPINK 

IF! PUTOUT! IOUT, 6) . EQ . 2 . 0 . AND. P .NE . 0 . 0 . AND. DISP .GT. D3 )  OUTPUT! 6) 

1'  CAUTION:  DISPLACEMENTS  AT  POINTS  GREATER  THAN  3D  AWAY  ', 

2*  FROM  THE  HQLE  MAY  OE  IN  ERROR' 

IF! PUTOUT! IOUT, 6) .EQ.2. )  WRITE! 6, 240) 

240  FORMAT! ///20X, ' DISPLACEMENTS' //'  DIST  ANGLE’ 10X,  '  U ' 10X,  'V'  ) 

IF( PUTOUT! IOUT, 6) .NE. 2 . )  GO  TO  270 
DO  260  J— 1 , NUMSTP 
DO  260  K— 1 , NUMPT 
ANGLE- (K-l ) ‘IAMG+ILOW 
DIST-! J-l ) ‘STPINK 

WRITE! 6, 250)  DIST, ANGLE , U (J, K ) ,V(J,K) 

250  FORMAT!  F7.3,F10.2,F13.6,F12.6) 

260  CONTINUE 
270  CONTINUE 
RETURN 
END 


•? 
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C 

SUBROUTINE  UNLODED( P, D1A, AI , BETA, STRESS, U, V ) 

CALCULATE  STRESS  DISTRIBUTION  AROUND  AN  UNLOADED  HOLE 

COMMON/THREE/  I ANG,  I  LOW ,  I  HIGH,  STPI NK,  NUMSTP 
REAL  I ANG, I LOW, I HIGH 

DIMENSION  STRESS ( 3,7,73),U(7,73),V(7,73),AI(3.3) 

DIMENSION  WORK( 5 ) ,COEF(5) ,RTR(4) , RTI(4) 

COMPLEX  Rl , R2 , COMPLX , XI 1 , XI 2 , COM1 , COM2 , DENI , DEN 2, PHI I ,  PHI2 
COMPLEX  Z,Z1,Z2,P1,P2,Q1,Q2 

CALCULATE  COMPLEX  PARAMETERS 

INITIALIZE  COMPLEX  NUMBER:  SQRT{-1.) 

COMPLX* ( 0. , 1 . ) 

NUMCO-4 

COEF(  I )  **AI  (2,2)  *  1000000 
COEF(  2)—  2.*AI(2, 3)  *1000000 
COEF( 3  )* (  2 . * AI ( 1, 2)+AI(3, 3) ) *1000000 
COEF( 4 ) *-2 . *AI (1,3) *1000000 
COEF  (  5  )  =*AI  ( 1 ,  1)  *1000000 
CALL  ROOTS  (COEF,  WO  RK.NUflCO,  RTR,  RTI,  IE) 

Rl*RTR(l )+COMPLX*RTI(l) 

IF(  RTI  (  2 )  •  GT  .0.0)  F  l=*RTR(  2  )+COMPLX*RTI  {  2  ) 

R2=8RTR(  3  )+COMPLX*RTI  (  3  ) 

IF( RTI ( 4 ) . GT .0.0)  R2»RTR( 4 )+COMPLX*RTI ( 4 ) 

P1-AIU, 1 ) *Rl *R1+AI (1,2) -AI (1,3)*R1 
P2-AI (1,1) *R2*R2+AI ( 1 , 2 )-AI ( 1 , 3 ) *R2 
Ql-AI ( 1 , 2 ) *R1+AI ( 2 , 2 ) /Rl-AI ( 2 , 3 ) 

Q2=*AI  (1,2  ) *R2+AI  (2,2  )  /R2-AI  (2,3) 


PI *3. 1415926535 
BETA=*BETA*  PI  / 1  BO .  0 

NUMPT=*(  ( I  HIGH- 1  LOW )  / 1  ANG )  + 1 

DC  20  JJ*1, NUMSTP 
DO  10  NtJ=l ,  NUMPT 

U(JJ,NN)-0.0 
V  ( JJ ,  NN )  =«0 . 0 
NNN*NN-l 

THETA8- ;  4fl*I  ANG+I LOW )  *PI/  1 30 . 0 
RADIUS=J  J  J*STP  I  NK+DIA/  2 . 0 

CALCULATE  X  £.  Y  COORDINATES  OF  POINTS  AROUND  UNLOADED  HOLE 

Xs8  RAD  I  US*  COS  (  THETA) 

Y=RADIUS*S IN (THETA) 

CALCULATE  LOCATION  PARAMETERS  FOR  UNLOADED  HOLE  EQUATIONS 

Z1=X+R1*Y 

Z2=X+R2*Y 

Z=X+COMPLX*Y 
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C 


MAPPING  FUNCTION 


XIl-CSQRT(Zl*Zl-DIA*DIA/4.-Rl*Rl*DIA*DIA/4.  ) 
XI2-CSQRT(Z2«Z2-DIA*DIA/4.-R2*R2*DIA*DIA/4. ) 

CHOOSE  CORRECT  SIGN  OF  CSQRT 

XI1-Z1/XI1 
XI2-Z2/XI.2 

IF(  REAL  (XI 1  )  .LT.-. 00001)  XI1»-1.*XI1 
IF  (  REAL(  XI 2  )  .LT.-.  00001)  XI 2— I.*XI2 

XI1*1 .-XII 
XI2-1 .-XI2 

CALCULATE  PHI  PRIME 

COM 1”R2*SIN( 2 . *BETA)+2 . *COS( BETA) *C0S ( BETA)+COMPLX* ( 2 . *R2* 
1  SI N ( BETA) *S IN( BETA ) +S IN( 2 . *  BETA ) ) 

C0M2“R1 *S IN( 2 . *  BETA )  >2 . *COS ( BETA) *COS( BETA) +COMPLX* ( 2 . *R1« 
1  SI  N  ( BETA ) *S IN( BETA )+S IN( 2 . *  BETA ) ) 

DENI  =«2  .  *OIA*  (  R1-R2  )  *  (  1  .  +COMPLX*Rl  ) 

DEN2=2 . *DIA* { R1-R2 ) * ( 1 . +COMPLX*R2 ) 

PH II—  COMPLX*P*DIA*COMl*XIl/{2.*DENl  ) 
PHI2=COMPLX*P*DIA*COM2*XI2/<2.*DEN2) 

CALCULATE  STRESSES  AROUNO  HOLE 

STRESS  (1.  JJ.  NN)=-P*COS(  BETA)  *COS(  BETA)  +2.  "REAL (R1*R1*PHI1+ 

1  R2*R2*PHI2) 

STRESS  (  2  #  JJ , NN)  '3'P*S  IN  (  BETA )  *S  I M(  BETA )  +  2  .  *  REAL  (  PHII  +  PHI2  ) 
STRESS(3.  JJ,NN)*P*S IN( BETA) *CCS( BETA) -2.* REAL (R1*PHI1  + 

1  R2*PHI2 ) 

CALCULATE  DISPLACEMENTS 

XI1=1 .-XII 
XI2-1 . -XI 2 

XI1=Z1/XI1 
XI 2  =  Z  2/Xl 2 

DF,N1  =  16  .  *  { R1-R2 ) * (  Zl+XI  I ) 

DEN2=16.* (RI-R2)*(Z2+XI2) 

PHI 1=-P*DIA*  DIA* ( COMPLX+R1 ) *COMl/DENl 
PHI2=P*DIA*DIA*  ( COMPLX+R2  )  *COM  2/DFN2 

U  ( J  J  ,  N N )  =  2  .  *  R EAL  (  PI  *  PHI  I +P  2*P!1 1 2  ) 

V( JJ,NN)=2 .  *  REAL  (  rJl  *  PH  I 1+Q2  *  PII 12  ) 

C 

10  CONTINUE 
20  CONTINUE 
C 

c 

RETURN 

END 
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SUBROUTINE  LOADED ( P,  DIA, S, ALPHA, STRESS, U,  V )  J 

CALCULATES  STRESS  DISTRIBUTION  AROUND  A  LOADED  HOLE  j 

ASSUMING  A  COSINE  BOLT  LOAD  DISTRIBUTION  | 

COMMON/ TWO/IOUT( 15) , NUMPLY, NUMMAT, ANG (8 )  , PLYTHK(O) ,MATID(3)  ] 

COMMON/THREE/ IANG, ILOW,IHIGH, STPINK , NUMSTP  j 

i 

REAL  IANG,  ILOW,  HUGH  1 

COMPLEX  R1 , R2 , COMPLX , Z , Z1 , 22 , CPOS ( 50), CNEG ( 50 ) , CZERO, CM ,  I 

1AKI , AK2 , XI 1 , XI2 , PHIL , PHI 2 , COMI , COM2 , XXI 1 , XXI 2  j 

COMPLEX  CHECK 1 , CHECK 2, PI,P2,Q1,Q2 


COMPLEX  Al (5 0 ) , A2 ( 50 ) 

DIMENSION  AMATRX (4,4), BMATRX (4,4) , STRESS( 3 , 7 , 73 ) 
DIMENSION  U( 7 , 7 3 ) ,V(7,73),S(3, 3) 

DIMENSION  WORK(5) , COEE ( 5 ) , RTR( 4 ) , RTI ( 4 ) 

INITIALIZE  COMPLEX  NUMBER:  SQRT(-1.) 

COMPLX* ( 0 . , 1 . ) 

CALCULATE  COMPLEX  PARAMETERS 
NUMCO=4 

COEF( 1)=S(2, 2)*1000000 
COEF(2)=-2.*S(2, 3) *1000000 
COEF(3)={2 .*S(1, 2)+S(3, 3) )*1000000 
COEF( 4 ) =-2 .  *S ( 1 , 3 ) *1000000 
COEF( 5 )=S( 1,1) *1000000 

CALL  ROOTS (COEF, WORK, NUMCO, RTR, RTI, IE) 

R1=RTR( 1 )+COMPLX*RTI ( 1 ) 

IF ( RTI ( 2 ) .GT.0.0)  R1=RTR( 2 )+COMPLX*RTI (  2  ) 

R2»RTR( 3 ) +COMPLX*RTI ( 3 ) 

IF ( RTI ( 4 ) . GT .0.0)  R2=RTR( 4 )+COMPLX*RTI (4) 

Pl=S(l, 1)*R1*R1+S(1, 2)-S(l, 3)*R1 
P2=S( 1, 1 )*R2*R2+S( 1 , 2)-S( 1,3) *R2 
Ql=S(l,2)*Rl+S(2,2)/Rl-S(2,3) 

Q2=S(1, 2)*R2+S(2,2)/R2-3{2, 3) 


PI=3.141592G535 

TIUCK=0. 

DO  10  N=l, NUMPLY 
THICK=THICK+PLYTHK (N ) 
10  CONTINUE 
P=4 . 0*P/PI 


A  COS I ME  LOAD  DISTRIBUTION  OVER  HALF  OF  HOLE  AT  AN  ANGLE 
ALPHA  TO  X  AXIS 

CALCULATE  COMPLEX  CONSTANTS 

PI2-PI/2.0 
M— 1 

20  CONTINUE 
M«=M+1 

IF(M.EQ.l)  GO  TO  40 
30  CONTINUE 

C1=SIN( ( M- 1 )  *  P 1 2 ) / ( 2  M  M- 1 ) ) 

C2=*SIN(  (M+l)*PI2)/(2*(fl+l)  ) 

C3=SIN( (M-I )*(-PI2) )/ (2* (M— 1 ) ) 

C4=SIN( (M+l ) * ( -PI2 ) )/(2*(M+l)) 

C5=*COS  (  (M-l  )*PI2  )  /  (  2*  (M-l )  ) 

C6=*COS  (  (M+I )  *PI2  )  /  (  2*  (M+l )  ) 

C7=*COS  (  (M-l  )*(-PI2)  ) /  (  2*  (M-l )  ) 

CB=»COS(  ( M+l  )*(-PI2)  )/  (2*  (M+l )  ) 

CM=»P*  (  ( C1+C2-C3-C4 )  -COMPLX*  ( -  C5-C6+C7+CB )  )  /  (2.0*PI) 
IF(M.EQ.O)  CZERO=CM 
iF(M.GT.l)  CPOS(M)=CM 
IF(M.LT.-l)  MN=-1*M 
IF(M.LT.-l)  CNEG(MH)=CH 
IF(M.LE.O)  GO  TO  50 
M— 1*M 
GO  TO  30 
40  CONTINUE 
C1=PI2 

C2=SIN(2.* (PI2) )/4. 

C3=*SIN (2.*(-PI2)  )/4. 

C4=SIN(PI2)*SIN(PI2)/2. 

C5=SIN(-Pl2)*SIW(-PI2)/2. 

CM=P*(  ( C1+C2-C3  )-M*COMPLX*  (  C4-C5 )  )/(2.*PI) 

IF(M.EQ.l)  CPOS(  1  )=*CM 
IF(M.EQ.-l)  CNEG ( 1 ) =CM 
IF(M.EQ.-1  )  GO  TO  50 
M=— 1 *M 
GO  TO  40 
50  CONTINUE 
M=I ABS ( M) 

IF(M.LT.40)  GO  TO  20 

TRANSFORM  COMPLEX  PARAMETERS  INTO  REAL  AND  IMAGINARY  PARTS 


S1  =  REAL(R1  ) 
S2=REAL(R2 ) 
T1=AIMAG ( R1 ) 
T2=AIMAG(R2 ) 


C  EQUATING  COEFFICIENTS  AND  SOLVING  FOR  CONSTANTS 

C 

DO  SO  M-1,45 
MN-M-1 

IF(MN.NE.O)  GO  TO  60 
DMATRX(l)«REAL(-CZERO*DIA/2. ) 

BMATRX(2)«AIMAG(-CZEKO*DIA/2. ) 

GO  TO  70 
60  CONTINUE 

BMATRX( 1 ) “REAL (- CPOS ( MN )*DIA/(2.*( MN+1 ) ) ) 

BflATRX(  2  ) “AIM AG  ( -CPOS{  MN )  *DIA/  (  2  .  *  ( MN+1 )  )  ) 

70  CONTINUE 
MN-M+1 
MNEG— 1*MIJ 

BMATRX*  3  )*REA7j(  -CNEG(MN)  *DIA/  (  2  .  *  (  MNEG+1 )  )  ) 
BMATRX(4)-AIMAG(-CNEG(MN)*DIA/  (2.*(MNEG+1  )  )  ) 
AMATRX(1,1)-T1+1. 

AMATRX( 1,2 )“S1 
AMATRX ( 1. 3 )*T 2+1. 

AMATRX ( 1 , 4)“S2 
AMATKX( 2 , 1 )»S1 
A’lATRX  l  2  ,  2  )“-Tl-l  • 

AMATRX { 2/3 )“S2 
AMATRX(2.4)»-T2-1  . 

AMATRX  (  3  ,  1 )  “1  .  -TI 
amatr:.(3, 2)— SI 

AMATRX ( 3 , 3 ) - I . -T2 

AMATRXl 3 , 4 )*-S2 

AMATRX ( 4 , 1 )-Sl 

AMATRX(  4 ,  2  )"*1  .  -T1 

AMATRX (4, 3) -S2 

AilATRX(  4 , 4 )“1  .  -T2 

CALL  SIMULT(AMATRX,B;1ATRX,4,J) 

IF( J . EQ . 1 ) OUTPUT ( 6 ) '  SIMULT  CALCULATES  A  SINGULAR  SET  OF  EQS 
A1  (M) =BMATRX( I ) +COMPLX* BMATRX ( 2 ) 

A2  ( M )  “15MATRX  (  3  )  >-COMPLX*nMATRX(  4  ) 

80  CONTINUE 
C 

PX“2  •  *  PI*  AIM  AG  ( COMPLX+  DIA*CNLG  (  1  )  /  2  .  ) 
PY-2.*PI*REAL(COMPLX*OIA*CNEG( I )/2. ) 

C 

AMATRX ( l ,  I  )=*T1 

AMATRX(1, 2)-Sl 

AMATRX (1,3) “T  2 

AMATRX ( 1 , 4 )“S2 

AMATRX ( 2 , 1 )“0 .0 

AMATRX(2,  2)  =  1  .0 

AMATRX (2,3) “0 • 0 

AMATRX (2, 4 )=l .0 

AMATRX(  3, 1 )“2  .  *'5L*T1 

AMATRX (3, 2  )=SI*S1-T1*T1 

AMATRX ( 3 ,  3  )”2 • *S2*T2 

AMATRX( 3 , 4 )“S2*S2-T2*T2 

AMATRX ( 4 , 1 )°-Tl/(Sl*51+Tl*Tl ) 

A1ATRX (4,2)=S1/(S1 *31 +T1 *T1  ) 

AMATRX(4, 3)*-T2/(S2*52+T2*r?) 

A/1ATRX ( 4 , 4  )  “S2/  ( S2*S2+T2*T2  ) 

BMATRX ( 1 )  =  PX/ ( 4.* PI) 

BMATRX ( 2 ) =-? Y/ ( 4 . *  P I ) 

BMATRX (3)“(S(1,2) *PY+3 (1,3)*PX)/(4.*PI*S(1,1)) 

BMATRX ( 4 )“-(S(l,2) "PX+S( 2,3) *PY)  / (4.*PI*S(2, 2) ) 

CALL  SIMULT( AMATRX, BMATRX, 4, J) 

IF( J. EQ. 1 ) OUT PUT (6 ) '  SIMULT  CALCULATES  A  SINGULAR  SET  OF  EQS 
C 

AK  l*n:iATRX  ( l )  +COMPLX* BMATRX  (  2  ) 

AK2-BMATRX { 3 ) +COMPLX*  Bt  JATRX ( 4 ) 
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C 


NUMPT-  (( I HIGH- I LOW ) / I ANG ) + 1 
ALPHA— ALPHA*  P I /  1  00 . 0 
ALPH— ALPHA 
DO  150  JJ-l.tlUMSTP 
DO  140  NN-l.WUMPT 

U(JJ,NN)=0.0 

V(JJ,NN)=0.0 

NNN=HN-1 

JJJ=JJ-1 

THETA=  (NUN* I ANG+I LOVJ )  *PI/ 180.0 
RADIUS=J JJ*3TPI NK+DIA/ 2 . 0 

CJVLCULATE  X  AND  Y  COORDINATES  OF  POINTS  AROUND  LOADED  HOLE 

X=RADIUS*COS ( THETA+ALPHA) 

Y=RAD I US * S I N ( T HETA+ALPHA ) 

CALCULATE  PARW1ETERS  FOR  LOADED  HOLE  EQUATIONS 

2 1=X+R1 *Y 
Z2=X+R2*Y 
Z=X+COMPLX*Y 

MAPPING  FUNCTION 


XXIl=CSQRT(Zl*Zl-DIA*DIA/4.-Rl*Rl*DIA*DIA/4. ) 
XXI 2=CSQRT( Z2*Z2-DIA*DIA/ 4 . -R2*R2*DIA*DIA/ 4 . ) 

CHOOSE  CORRECT  SIGN  OF  CSQRT 

90  CONTINUE 
XI1=Z1+XXI1 
XI2=Z2+XXI2 

XI 1=2 . *XI1/ ( DIA* ( 1 . -COMPLX*R^ ) ) 

XI  2=2 .*XI2/ ( DIA* ( 1 .-C0MPLX*R2) ) 

C0X1=REAL (XII) *REAL (XII) + AIM AG (XII) * AIM AG (XII) 
COX2=REAL(XI2 ) *REAL(XI2 )+AIMAG(XI2 )*AIMAG(XI2) 
IF(COXl.GE. .99999)  GO  TO  100 
XXI1—XXI1 
GO  TO  90 
100  CONTINUE 

IF(COX2.GE. .99999)  GO  TO  110 
XXI2—XXI2 
GO  TO  90 
110  CONTINUE 
XXI1=XI1 
XXI 2=XI 2 

CJ1LCULATE  PHI  PRIME 


C0M1= ( 0 . ,0.  ) 

COM2=(0 . ,0. ) 

DO  120  M=1 , 45 

C0M1=C0M1+M*A1 (M)*XI1**(-1*M) 
COM2=COM2-K1*A2(M) *XI2** (-1*M) 
120  CONTINUE 
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CHECK  SIGN  OF  CSQRT 

XIl=CSQRT(Zl*Zl-DIA*DIA/4.-DIA*DIA*Rl*Rl/4.  ) 

XI2=CSQRT(  Z2*Z2-0'IA*DIA/4  . -DIA*  DIA*  R2*R2/4 .  ) 

CHECK1=Z 1 /XI 1 
CHECK2-Z2/XI2 

IF ( REAL( CHECK1 ) .LT.-. 00001)  XI1=-1 .*XI1 
IF(REAL(CHKCK2) .LT.-. 00001 )  XI2=-1 .*XI2 
PH 11= ( AKI-COM1 ) /XI 1 
PHI2= ( AK2-COM2 ) /XI2 

CALCULATE  STRESS  COMPONENTS  IN  LAMINATE  AT  COORDINATES  X.Y 

STRX=2  .*REAL(R1*R1*PHI1  +  R2*R2*PIII2) 

STRY=2 .* REAL (PHI 1+PHI2) 

STRXY=-2 . *REAL( R1*PHI1+R2*PHI2 ) 

STRESS ( 1 , JJ,NH)=STRX*COS(  ALPH) *COS ( ALPU ) +STRY*S IN( ALPH ) * 

1  SIN(ALPH)-2.*STRXL’’*SIN(  ALPH)  *COS( ALPH) 

STRESS(2, JJ,NN) =STRX*SIN( ALPH) *SIN( ALPU )+STRY*C05( ALPH)* 

1  COS ( ALPH )+2. *STRXY*SIN( ALPU )*COS( ALPH) 

STRESS ( 3 , JJ , NN)=STRX*SIN( ALPH) *COS( ALPH ) -STRY*SIN( ALPH ) * 

1  COS ( ALPH )+STRXY* (COS v ALPH ) *COS(  ALPH  )- 

2  SIN(AI.PH)  *SIN(ALPH)  ) 

CALCULATE  DISPLACEMENTS 

XI1=XXI1 
XI 2=XXI 2 
C0M1= ( 0 . ,0. ) 

COM  2= ( 0 . ,0. ) 

DO  130  M= 1 ,45 

COMl=COMl+Al (M) *XI1** ( - 1 *  M ) 

COM2=COM2+A2(M) *XI2** (-1*M) 

130  CONTINUE 

XXI l=CLOG (XII  ) 

XXT.2=OLOG(XI2  ) 

PHI1=AK1*XXI 1+COM1 
PHI2=AK2*XXI2+COM2 
U( JJ,UN)=2 .*REAL(P1*PHI1+P2*PHI2) 

V(JJ,NN)=2 . *REAL( Q1*PHI1+Q2*PI1I2 ) 

140  CONTINUE 
150  CONTINUE 

RETURN 

END 
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SUBROUTINE  PLYSTR( IFAIL) 

TRANSFORMS  LAMINATE  STRAINS  TO  PLY  STRESSES/ STRAINS  BY  ASSUMING 
CONSTANT  STRAIN  THROUGH  THE  THICKNESS 

COMMON/ ONE/ El ( 3 ) , E2 ( 3 ) , G12 ( 3 ) , VI2 ( 3 ) 

COMMON/TWO/ I OUT( 15) ,NUMPLY, NUMMAT, ANG{8) ,PLYTHK(8) ,MATID(8) 
COMMON/' THREE/ I ANG,  I  LOW ,  IHIGH ,  STPINK,  NUMSTP 
COMMON/ EIGHT/ STRESS ( 3,7,73), STRAIN ( 3,7,73) 

COMMON/NINE/ STR1 (8 , 7,73), STR2 (8,7,73), STR1 2(8,7,73) 

REAL  IANG, I LOW , IHIGH 


i 


* 


STRAINS  PER  PLY 
MOVE-O 

NUMPT-( ( IHIGH-ILOW) /IANG)+1 

IF( PUTOUT( I OUT, 7 ) , EQ . 2 . )  WRITE(6, 10) 

10  FORMAT(///20X, ‘STRAINS  PER  PLY*,//'  DIST  ANGLE  PLY', 

1*  STRAIN  1  STRAIN  2  SHEAR  STRAIN*/) 

20  CONTINUE 

DO  40  JJ-1, NUMSTP 
DO  40  NN— 1 , NUMPT 
DO  40  L-l.NUMPLY 
D«ANG(L)*3. 1415926535/180.0 
STRANX-STRAIN(1, JJ,NN) 

STRANY-STRAIN( 2 , J J , NN ) 

GAMA-STRAI N ( 3 , J J , NN ) 

STRANl-STRANX*COS(D)*COS(D) 

STRAN2-STRANY*  SIN(D)*SIN(D) 

GAMA1 2 “GAMA* SIN(D)*COS(D) 

STR1(L, JJ,NN)-STRAN1+STRAN2+GAMA12 
STRAN1-STRANX*SIN(D)*SIN(D) 

STRAN2-STRANY*COS ( D) *COS( D ) 

GAMA12— 1  .*GAMA*SIN(D)*COS(D) 

STR2(L, JJ,NN)-STRAN1+STRAN2+GAMA12 
STRAN1— 2  .  *STRANX*SIN(  D )  *COS(  D ) 

STRAN2-2 .*STRANY*SIN(D)*COS(D) 

UAMA1 2-GAMA*COS ( D ) *COS ( D ) -GAMA*  SIN(D)*SIN(D) 

STR12(L, JJ,NN)«STRAN1+STRAN2+GAMA12 

ANGLE- (NN-1 }*IANG+IL0W 
DIST-( JJ-1 ) *STPINK 

IF( PUTOUT( IOUT, 7 ) ,EQ . 2 . )  WRITE( 6 , 30)DIST, ANGLE, ANG( L) , 
1STR1(L, JJ.NN) ,STR2(L, JJ.NN) ,STR12(L, JJ.NN) 

30  FORMAT( F6 . 3 , 2F8. 2 , 3F1 2 .6) 
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STRESSES  PER  PLY 

IF( PUTOUT( I OUT, 0) .EQ.2. )  WRITE(6, 50) 

50  F0RMAT(///20X, 'STRESSES  PER  PLY',//'  DIST  ANGLE  PLY', 

1'  STRESS  1  STRESS  2  SHEAR  STRESS'/) 

DO  70  JJ-l.NUMSTP 
DO  70  NN-l.NUMPT 
DO  70  L-l.NUMPLY 
M-MATID(L) 

V21-V12  (M)*E2(M)/El(M) 

DEN-1. -V12(M)*V21 
ABC-STR2 ( L, JJ , MN) /DEN 
BCA-STRKL,  JJ.NN) 

STR1(L, JJ,NN)»E1(M)*STR1(L, JJ, UN) /DEN+V12 (M) *E2 ( M) *ABC 
STR2(L,  JJ.NN) -VI 2  (.'‘1)  *E2(M)  *BCA/DEN+E2 (M) *ANC 
STR12 ( L, JJ, NN) -STR12 ( L,  JJ, MU) *G1 2 ( M) 

ANGLE- (  Ntl-l  )  *1 ANG+ILOW 
DIST— ( JJ-1 ) *STPINK 

IF(  PUTOUT(IOUT,8)  .EQ.2.)  WRITE<  6,  60)  DIST,  ANGLE, ANG(L) , 
1STR1(L,  JJ.NN) ,STR2(L, JJ.NN) ,STR12(L, JJ.NN) 

60  FORMAT( F6 . 3 , 2F8 . 2 , 3F12 . 2 ) 

70  CONTINUE 
MOVE-1 

IF( IFAIL. EQ . 1 )  GO  TO  20 
80  CONTINUE 

RETURN 

END 


SUBROUTINE  FAILURE 

POINT  STRESS/STRAIN  ANALYSIS  FOR  FAILURE  USING  UNIDIRECTIONAL 
MATERIAL  ALLOWABLES 

COMMON/ TWO/ I OUT ( 15), NUMPLY, NUMMAT, ANG( 8 ) , PLYTHK( 3 ) , MATID( 8 ) 
COMMON/THREE/IANG, ILOW, THIGH, STPINK, NUMSTP 
COMMON/ FOUR/ PX , PY , PXY , P , PW , ALPHA , BETA , DI A , CORRECT 
COMMON/ FIVE/ FXT ( 3 ) , FXC( 3 ) , FYT ( 3 ) , FYC ( 3 ) , FXY( 3 ) , IFAIL 
COMMON/NX NE/ STR1 (8,7,73),STR2(8,7„73), STR1 2(8,7,73) 

REAL  IANG, ILOW, I HIGH 

DIMENSION  PLYFAII.(  3,8),  FAI LS  (  3 , 0 ) ,  RATIO(  3,8),  PLYRATO(  3,8) 

IF(PUTOUT(IOUT,9) .NE.2.)  GO  TO  40 
IF( IFAIL.GT . 2 )  GO  TO  20 
WRITE( 6,10) 

10  FORMAT(///2 OX, ‘FAILURE  CRITERIA  PER  PLY',//'  DIST  ANGLE  ', 
1'  PLY' , 12X,  '  FAILURE  NUMBERS' /35X, ' 1  2  SHEAR') 

GO  TO  40 
20  CONTINUE 
WRITE( 6 , 30) 

30  FORMAT(///20X, 'FAI LUBE  CRITERIA  PER  PLY',//'  DIST  ANGLE', 
1'  PLY  FAILURE  FAILURE  RATIOS ’/ 3 2X ' NUMBER 

2*  1  2  SHEAR’) 

40  CONTINUE 
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CHECK-0 . 

KKX-1 

F2-0.0 

F3-0.0 

NUMPT-( ( IHIOH-ILOW) /IAN  )  +  l 
DO  170  JJ— 1 , NUMSTP 
DO  170  KK-l.NUMPT 
SIQ-1.0 

DO  160  I I— 1 , NUMPLY 

X-STR1 ( 1 1 ,  J  J , KK) 

Y-STR2 ( 1 1 ,  J  J , KK) 

XY-STF.12(  II,  JJ,KK) 

MATH— HAT  I  D(  II) 


GO  TO  (50,50,60,70,80)  1  FAIL 
50  CONTINUE 

MAXIMUM  STRESS/ STRAIN 
FX-FXT ( MATI I ) 

IF(X.LT.O.O)  FX-FXC(MATII) 
FY-FYT ( MAI  II) 

IF(Y.LT.O.O)  FY— FYC(MATII) 

Fl-X/FX 

F2-Y/FY 

F  3— XY / FXY ( MATI I ) 

GO  TO  90 


60  CONTINUE 

TSAI -HILL  CRITERIA 
FX-FXT (MATI I) 

IF(X.LT.O.O)  FX— FXC ( MAT II) 

FY-FYT (MATI I) 

IF(Y.LT.O.O)  FY— FYC ( MAT I I ) 

F1-X*X/ ( FX"FX)+Y*Y/ (FY*FY)-X*Y/ ( FX*FX)+ 
1  XY*XY/ ( FXY( MATII ) *FXY( MATI I ) ) 

RATIOX— (X/FX)/SQRT(F1 ) 
RATIOY«(Y/FY)/SQRT(Fl ) 

RATIOXY— ( XY/ FXY( MATI I ) )/SQRT(FJ ) 

GO  TO  90 


70  CONTINUE 

MODIFIED  TSAI-V/U  CRITERIA 

Fl-1 . / FXT( MATI I ) -1 . /  FXC ( MATII ) 

F2-1 ./FYT( MATII )-l ./FYC (MATI I) 

FI 1—1 . / ( FXT ( MATII ) *FXC ( MATII ) ) 

F22-1 ./ (FYT(MATII)*FYC(MATII) ) 

F66— 1 . / { FXY ( MATII ) * FXY ( MATII ) ) 
F1-F1*X+F2*Y+F1 1*X*X+F2  2*Y*Y+F66  ■>XY*XY 
FX-FXT (MATI I) 

IF(X.LT.O.O)  FX-FXC ( MAT II) 

FY-FYT (MATI I) 

IF(Y.LT.O.O)  FY-FYC(MATin 
FI-AQS(Fl) 

RATIOX- ( X/FX ) /SQRT( FI ) 

RATIOY- { Y/ FY ) / SQRT( FI ) 

RATIOXY- (XY/ FXY (MATI I) )/SQRT(Fl ) 

GO  TO  90 
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80  CONTINUE 


HOFFMAN  FAILURE  CRITERIA 

Fl-1  •  /FXT (  MATH )  -i  .  / FXC ( MATH  ) 

F2-1 .  /FYT(  MATH  )-l .  /FYC(  MATH ) 

Fll«l./  (  FXT  (MATII)*FXC(  MATH)  ) 

F22-1  .  /  (  FYT( MATII )  *FYC( MATI I )  ) 

F66-1 . / ( FXY( MATI t ) *FXY( MATI I ) ) 

F12— 1 ./(FXT(MATll)*FXC(MATII) ) 

F1«F1*X+F2*Y+F11*X*X+F22*Y*Y+F12*X*Y+F66*XY*XY 

FX-FXT(MATII) 

IF(X.LT,0.0)  FX-FXC(MATII) 

FY-FYT(MATII) 

IF(Y.LT.O.O)  FY-FYC(MATII) 

Fl-ABS(Fl) 

RATIOX-(X/FX)/SQRT(Fl ) 

RATIOY-( Y/FY)/SQRT(F1) 

RATIOXY-( XY/ FXY ( MATI I ) ) /SQRT{ FI ) 

GO  TO  90 


90  CONTINUE 

ANGLE- (KK-1 ) *IANG+ILOW 

DIST— ( JJ-I ) *STPINK 

IF( IFAI L. GT . 2 )  GO  TO  110 

I F ( PUTOUT ( I OUT ,  9 ) . EQ . 2 . )  WRITE(6.100)  DIST, ANGLE, ANG( II), 
1F1 1 F2»  F3 

100  FORMAT ( F7.3,2F10.2,3F10.3) 

GO  TO  130 
110  CONTINUE 

I F( PUTOUT ( I OUT, 9 ) . EQ . 2 . )  WRITE(6,120)  DIST, ANGLE, ANG( I I ) , 
1F1, RATIOX, RATIOY, RATIOXY 
120  FORMAT ( F7.3,2F10.2,4F10.3) 

130  CONTINUE 

AUTOMATIC  SEARCH  FOR  FAILURE  ROUTINE 

IF( S iG . EQ . 2 . )  FAILS(1,II)-F1 

IF( SIG . EQ •  2  •  )  FAILS ( 2, II}— F2 

IF(SIG.EQ.2. )  FAI LS ( 3  » 1 1 ) - F  3 

IF(SIG.EQ. 2. )  RATIO( III )— RATIOX 

IF( SIG. EQ. 2 . )  RATIO( 2,11) -RATIOY 

IF(  SIG.  EQ  •  2  .  )  RATIO(  3,11 )  —  RATIOXY 

IF( JJ . NE. 2 )  GO  TO  150 

PLYt’AI  L(  1,11 )— FI 

PLYFAIH  2,  II)-F2 

Pl,YFAIL(  3,  II)r»F3 

PLYRATO( 1 , I I )-RATIOX 

PLYRATO( 2 , I I ) -RATIOY 

PLYRATO( 3,11 ) —RATIOXY 

CIIK— CHECK 

IF( ABS( CHECK ) .  LT. ABS(Fl) )  CHECK-Fl 

I F( AOS ( CHECK ) . LT. AUS( F2 ) )  CHECK-F2 

IF (AOS (CHECK) .LT. AnS( F3) )  CHECK-F3 

IF{ CHECK . EQ. CHK )  GO  TO  150 

KKK-KK 

III-II 

SIG-2.0 

DO  140  -1-1,111 

DO  140  tl-1 ,  3 

FAI LS ( N, M)=PLYFAIL(N,M) 

RATIO (N, M) »PLYRATO( N, M) 

140  CONTINUE 
150  CONTINUE 
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160  CONTINUE 
170  CONTINUE 


IF ( CHECK. EQ. 0.0)  GO  TO  260 

IF(  IFAIL.  EQ.  1)  CORRECT- 1 . 0/ABS ( CHECK ) 

IF ( IFAIL. EQ. 2 )  CORRECT-1 . 0/ABS( CHECK ) 

IF( IFAIL. EQ. 3 )  CORRECT-1 . 0/SQRT( CHECK ) 

IF( IFAIL. EQ . 4 )  CORRECT-1 . 0/SQRT( CHECK ) 

IF( IFAIL. EQ. 5 )  CORRECT-1 . 0/SQRT( CHECK) 

IF( PUTOUT( I OUT, 10 ) . NE. 2 . )  GO  TO  250 
IF (CORRECT.LT. .990. OR. CORRECT. GT. 1 .001)  GO  TO  250 
ANGLE- (KKK-1 )*IANG+ILOW 

IF( IFAIL.GT. 2)  GO  TO  210 
WRITE ( 6 , 100 )  PX,PY,PXY,P 

mo  FORMAT(///20X. 'AUTOMATIC  SEARCH  FOR  FAILURE* // 2 5X, 

1‘ FAILURE  STRESSES '/1SX, ' PX ' , 10X , ' PY ' , 1 OX , ' PXY' , 10X, ' P* / , 
211X.4F12.2//9X, 'DIST  ANGLE  PL*  >  FAILURE  NUMBERS' / , 

336X, ' 1  2  SHEAR') 

DO  200  I-i , NUMPLY 

WRITE ( 6 , 190 )  STPI NK, ANGLE , ANG ( I  )  , FAI LS ( 1 , 1 ) , FAI LS ( 2 . I ) . FAILS'  3.1) 
190  FORMAT ( FI l . 3 , F9 . 2 , F10 . 2, 3F0 . 3 ) 

200  CONTINUE 
GO  TO  260 

210  CONTINUE 

'WRITE ( 6 , 220  )  PX.PY.PXY.P 

220  FORHAT(///20X, 'AUTOMATIC  SEARCH  FOR  FAILURE '// 25X, 

1 'FAILURE  STRESSES' /1SX, ' PX’ , 10X, '  PY'  , 10X, ' PXY’ , 10X, ' P' / , 
211X.4F12.2//5X, ’  DIST  ANGLE  PLY  FAILURE 

39X, ' FAILURE  RATIOS' / 3 5X,  'NUMBER' . 7X, ' 1 ' , 9X , *  2 '  , 7X,  ' SHEAR*  ) 

DO  740  1-1, NUMPLY 

WRITE( 6,2  30)  STPINK, ANGLE,  ANG( I ) ,FAILS(1,  I )  ,  RATIOl  1,1)  , 

1  RATIO( 2,1), RATIO( 3,1) 

230  FORMAT ( Fll . 3 , F9 . 2, F10. 2 , 4F10 . 3) 

240  CONTINUE 
250  CONTINUE 


260  CONTINUE 

IF( PUTOUT ( IOUT ,10) .ME. 2 . )  CORRECT-1 .0 

RETURN 

END 
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sun rckjt i ne  s imult (a.u.u.xs) 

test  tor  algorithmic  a::m\Rm  added  01/10/T9 

DIMENSION  ACL).U(L) 

MAC  Hi;  MU  EPSILON  FOR  CYDER  SINGLE  PRECISION 
DATA  GPS/ 7 • 11E-15  f 
TOL-3.IG*r.PS»(M-l ) 

3ETA-0 . O 

,JJ— M 

01  r»5  j-i,n 
JY*J+l 

IJT'5A*0.0 
IT-JJ-J 
DO  30  I»J,U 
IJ-IT+I 

IF(M$S(  ni'SA )- A0S(  A(  IJ)  )  )  20,  30,  30 
20  JIGA-M  IT) 

IMAX* I 

30  CONTINUE 

IF(  A0S(  RIGA )  .GT.  OETA)  3ETA«A3S(  QIGA) 

IFlAHS( niGA)-TOL*UETA)  35,35,40 
35  XS-I 
RETURN 

40  IW+.N*(J-2) 

IT-IMAX-.J 
L'O  50  K*J,N 

11- Il+N 

12- I1+IT 
SAVE-A( I 1 ) 

A(I1)-A(I2) 

A( 12) -SAVE 

50  \(  I  l)-A(  ID/BIGA 
3AVE-3(  I  MAX) 

U(I  lAX)-0( J) 

B( J)-S\VK/BIGA 
IF(J-M)  55,70,55 
55  INS-N*(J-1) 

DO  65  IX-.TY,H 

IXJ-IQS+IX 

IT-J-IX 

DO  60  JX-JY.N 

IXJX-U*(JX-l)+IX 

J JX-IXJX+IT 

60  A(  IXJX)-A(  IX.7X ) - ( A(  IXJ)*A(JJX)  ) 

65  o( ix)»n(ix)-(n(J)*A(ixj) ) 

70  NY-N-1 
IT»N*N 

DO  00  J-l.UY 
IA-IT-J 
1 3*  N-J 
IC»U 

DO  BO  K*1 , J 

B(  IB)»lt(ID)--A(  IA)  *D(  IC) 

IA-IA-M 
80  IC-1C-I 
RETURN 
END 
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SUBROUTINE  ROOTS ( XCOF , COF , H, ROOT R,  ROOT  I,  IF.R) 
DIMENSION'  XCOF(M)  ,COF(M)  ,ROOTR(M)  ,  ROOTl(M) 

DOUULE  PRECISION  XO , YO, X, Y, XPR , YPR , UX , UY, V , YT , XT , U 

1  , XT2 , YT2 , SUMSQ, DX , DY, TEMP,  ALPHA,  FI 

2  ,  RMPRF.C ,  TOL 

C  RELATIVE  MACHINE  PRECISION  (TEST  FOR  'ALMOST  ZERO') 

DATA  RMPREC/I .00-14/  .TOL/l.OD-4/ 

IFIT-0 

N-M 

IER-0 

IF( XCOF ( N+l ) )  10,30,10 
10  IF ( N )  50,20,50 
20  IER-1 

GO  TO  360 
30  ISR-4 

GO  TO  360 
40  IER-2 

GO  TO  360 

50  IF( N-36 )  60,60,40 
60  MX-H 
NXX-M+1 
N2-1 

KJ1  »  N41 
DO  70  L-l.KJl 
MT-KJl-L+1 
70  COF (MT) -XCOF (L) 

90  XO». 005001 01D0 
Y0=0. 0100010100 
IN-0 
90  X-XO 

XO— IO.DO*YO 
YO=»  — 1 0  .  D0*X 
X-XO 
Y*YO 
IN- I N+l 
GO  TO  no 
100  IFIT=*1 
XPR~X 
YPR=>Y 
110  ICT-0 
120  UX-O.DO 
UY-0 . 00 
V  =*0  .  DO 
YT-0.00 
XT-1  .00 
U=COF( N+l ) 

IF(D,\nS(U)  .LE.  RMPRF.C) GO  TO  200 
130  00  140  1=1, V 
L  -N-I+l 
TtMP-COF(L) 

XT2— X*XT-Y*YT 
YT2-X*YT+Y*XT 
U-U+TEMP+XT2 
V»V+TEMP*YT2 
FI -I 

UX-UX+FI+XT*TEMP 
UY=UY-FI*YT*TF,MP 
XT=XT2 
140  YT-YT2 

SUMSO=UX*UX+UY*UY 
1F(  SHMSQ.LF. .  RMPRl’C)GO  TO  2  30 
150  DX- ( V*UY-U*UX )  / SUMSQ 
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X*X+DX 

DY—  ( U*UY+V#UX  )  / SUMSQ 
Y-Y+DY 

160  IF(  DABS(DY)+  DABS(DX) . LE.TOL)  GO  TO  210 
170  ICT-ICT+1 

IF ( ICT-500 )  120,190,130 
130  IF(IFIT)  210,190,210 
190  IF( IN— 5 )  90,200,200 
200  IER-3 

GO  TO  360 

210  DO  220  L-l.NXX 
MT-KJ1-L+1 
RTEMP-XCOF(MT) 

XCOF ( MT ) *COF ( L ) 

220  COF(L.)-RTEMP 
ITEMP-N 

:j-nx 

NX* I TEMP 

IF(IFIT)  250,100,250 
230  IF(IFIT)  740,90,240 
240  X*XPR 
Y»YPR 
GO  TO  210 
250  IFIT-0 

260  IF ( DAB3 ( Y )- 1 . 0D-4*DABS ( X ) )  290,270,270 
270  ALP'IA-X+X 

SUMSQ-X*X+Y*Y 
N-N-2 
GO  TO  300 
230  X*0 . DO 
NX-MX-1 
NXX-NXX-1 
290  Y-O.DO 

SUMSQ-O.DO 

ALPUA*X 

N-N-l 

300  COF( 2 ) =COF ( 2 )  +AL.PH  A*COF ( 1 ) 

310  DO  320  L*2,N 

320  COF( L+l )*COF( L+l ) +ALPHA*COF( L)-SUMSQ*COF( L-l ) 
330  ROOTI (N2 )*Y 
ROOTR(N2 )*X 
N2-M2+1 

IF(SUMSQ.I  E. RMPREC)  GO  TO  350 
340  Y*-Y 

SUMSQ*0.D0 
GO  TO  330 

350  IF( N . GT . 0 ) GO  TO  90 
360  RETURN 
END 
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